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Abstract 

A  constitutive  equation  is  proposed  with  a  view  to  describing  the  rate 
dependent  mechanical  response  of  metals  at  high  temperatures.  The  equation 
is  of  the  endochronic  type  and  derives  its  physical  foundations  from 
deformation  kinetics.  Of  importance  is  the  fact  that  hardening  is 
associated  with  a  change  in  the  energy  barriers  brought  about  by  the 
inelastic  deformation  of  a  metal.  The  equation  is  used  to  describe  the 
results,  by  Ohno  and  his  associates,  of  experiments  on  the  creep  response 
of  metals  to  piece-wise  constant  stress  histories.  The  metal  in  the 
present  case  is  304  stainless  steel  at  600°  C.  It  is  shown  that  the  theory 
gives  analytical  results  that  are  in  agreement  with  exDeriment. 

Of  consequence  is  the  fact  that  the  constitutive  equation  applies  to 
three-dimensional  stress  or  strain  histories  and  is  thus  not  limited  to 

those  stress  histories  associated  with  creep. 

The  theory  is  also  extended  to  large  deformations.  This 
is  done  by  using  the  internal  variable  theory.  The  resulting 
constitutive  equation  is  a  statement  to  the  fact  that  the  Cauchv 
stress  is  a  quadratic  functional  of  the  relative  Finger  Tensor 
in  terms  of  an  intrinsic  time  wich  is  defined  in  the  text. 


1.  Introduction 


In  the  present  paper  we  develop  an  endochronic  theory  of  viscoplasticity 
which  accounts  for  the  history  of  strain  and  strain  rate  on  the  stress  response 
of  metals  at  high  temperatures.  The  theory  is  based  on  the  concepts  of 
endochronic  plasticity  (see  typically  Ref.'s  [1],  [2],  and  [3]},  however,  the 
increment  of  intrinsic  time  scale  is  no  longer  proportional  to  the  plastic 
strain  path  but  depends  also  on  the  rate  at  which  the  path  is  traversed.  The 
development  of  the  theory  is  dealt  with  at  length  in  the  subsequent  sections. 

The  resulting  constitutive  equation  is  used  to  analyze  the  creep  response 
of  304  stainless  steel  to  piece-wise  constant  shear  stress  histroies  at  600°C 
and  to  compare  the  results  with  the  experimentally  determined  creep  response  of 
the  same  material  at  this  temperature  as  reported  by  Ohno  et  als.  in  Ref  [4], 

The  constitutive  behavior  of  metals  at  high  temperature  where  the  strain 
rate  sensitivity  of  the  mechanical  response  cannot  be  ignored,  has  been  the  sub¬ 
ject  of  extensive  research  in  recent  years.  We  do  not  wish  to  give  in  this 
paper  an  exhaustive  review  of  the  literature  on  this  subject  but  merely  cite 
references  which  are  typical  of  the  enormous  amount  of  work  which  is  being  done 
ir  this  field.  In  this  context,  the  works  of  Chaboche  [5],  Krieg  [6],  Malvern 
[7],  Haisler  [8],  Bradley  [9],  Leckie  [10],  Krempl  [11],  Walker  [12],  Miller 
[13],  and  Hart  [14],  among  others  are  mentioned.  We  also  wish  to  cite  the  works 
of  Ohashi  et  als.  [15],  and  Murakami  and  Ohno  [16],  and  Ohno  et  als.  [4]  who 
have  been  enjoying  a  measure  of  success  in  using  a  creep  hardening  surface 
theory  in  describing  the  creep  response  of  metals  to  piece-wise  constant  stress 
histories,  a  subject  with  which  we  will  be  dealing  in  this  paper.  This  aspect 
of  the  mechanical  response  of  metals  has  given  rise  to  greater  difficulties 
than  say,  the  stress  response  to  piece-wise  constant  strain  rate  histories. 
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The  work  of  Krausz  [17]  and  his  co-workers  also  occupies  a  significant  position 
in  the  literature  because  of  its  more  fundamental  nature  in  that  a  microscopic 
theory  of  deformation  kinetics  is  used  to  gain  understanding  of  the  mechanical 
response  of  metals  at  the  bulk  level. 

In  our  initial  approach  to  the  subject,  with  specific  attention  to  the 
strain  (creep)  response  to  piece-wise  constant,  pure  shear  stress  histories,  we 
used  a  strictly  phenomenological  approach,  in  the  context  of  an  endochronic 
theory.  Specifically  in  one  dimension  we  probed  the  data  with  the  constitutive 
equation 


eP  =  /  J(z-z')  Jfr  dz' 


0- 


(1.1) 


Henceforth  convolution  integrals  such  as  the  one  appearing  on  the  right  hand 
side  of  eq.  (1.1)  will  be  given  in  symbolic  form  according  to  eq.  (1.1a): 


z  d  def 

J(z  -  z ' )  -j—  i  dz'  =  J(z)*ds 
Cl¬ 
in  equation  (1.1)  eP  has  the  same  connotation  as  eC, 


shear  strain  where 


(1.1a) 

i.e.,  it  is  the  inelastic 


•p  *c 
e  =  e  =  e 


(1.2) 


and  Uq  is  the  elastic  shear  modulus.  We  caution,  however,  that  other  defini¬ 
tions  of  ec  have  been  used  in  the  literature.  The  intrinsic  time  z  was  defined 
by  equation  (1.3)  where 


dz 


dc 

f(c,c) 


(1.3) 


•  •  £ 

and  c  =  |e  |  is  the  usual  fashion.  The  dependence  of  f  on  c  lends  "strain  rate 
sensitivity"  to  the  equation,  which  otherwise  would  be  strain  rate  insensitive. 


While  equations  of  the  type  (1.1),  (1.2)  and  (1.3)  were  shown  to  give  satis¬ 
factory  results  in  the  case  of  constant  strian  rate  histories  as  demonstrated  by 
Wu  and  Yip  [18,19]  and  Lin  and  Wu  [20,21],  we  found  that  these  equations  did  not 
prove  satisfactory  in  the  case  of  piece-wise  constant  stress  histories.  With 
specific  reference  to  the  data  of  Ref.  [4],  it  was  found  that  f  had  to  be  essen¬ 
tially  independent  of  z,  to  account  for  the  periodic  creep  response  to  piece- 
wise  constant  cyclic  histories.  Thus,  limiting  f  to  a  dependence  on  c  only 
resulted  in  a  gross  overestimate  of  the  creep  strain  under  cyclic  conditions. 
Correcting  f  so  as  to  match  the  data  gave  rise  to  oscillations  in  f  which  could 
not  be  accounted  for  by  means  of  equation  (1.3).  Furthermore,  the 
phenomenological  approach  did  not  give  any  hint  as  to  the  physical  mechanism(s) 
responsible  for  such  fluctuations  in  f.  To  overcome  the  difficulties  we  appealed 
to  the  theory  of  deformation  kinetics  in  the  context  of  the  internal  variable 
theory.  The  latter  is  treated  briefly  in  Section  2  while  the  former  is  repre¬ 
sented  in  detail  in  Section  3. 
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2 .  Internal  Variable  Theory 


Irreversible  thermodynamics  of  internal  variables  is  now  a  well  established 
field  so  we  proceed  to  give  a  very  breif  outline  of  the  theory  for  the  sake  of 
completeness.  We  limit  ourselves  to  small  strain  fields.  In  the  Helmholtz  for¬ 
mulation  the  thermodynamic  state  is  described  by  the  free  energy  density  i|i  which 
is  a  function  of  the  strain  tensor  6  the  temperature  T  and  n  interval  variables 
qr  which,  for  thermomechanical  processes,  are  tensors  of  the  second  order. 

M 

The  stress  6  and  the  entrpy  density  ^  are  then  given  by  equations  (2.1)  and 

(2.2) 


(2.1) 


3(f) 

V  -  ST 


(2.2) 


In  the  case  of  a  spatially  uniform  thermal  field  the  condition  of  positive  rate  of 
irreversible  entropy  gives  rise  to  the  inequality 


3  ip  •  r  •  r 

— -  •  q  >  0  ,  i  q  i  *  0 

3q 


(2.3) 


r,  not  summed. 

In  the  linear  version  of  the  endochronic  theory  i|>  is  a  quadratic  function  of 

.  .  r 

its  arguments  and  the  evolution  equations  for  q  have  the  form 


3^ 

3q 


dqr 


_  +  b 

r  «r 


dz 


=  0 


(2.4) 


r  not  summed,  where  n  such  equations  exist,  one  for  each  variable  q  ,  and  b 


are  positive  definite  tensors  of  the  fourth  order. 
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When  isothermal  conditions  prevail  and  the  solid  is  initially  isotropic  and 

the  tensors  b  are  constant,  then  equations  (2.1)  and  (2.4)  combive  to  give  rise 
r 

M 

the  integral  constitutive  equations 

s  =  2u(z)*de  (2.5) 


o  =  3K(z)*de 


(2.6) 


where  s  is  the  stress  deviator,  e  the  strain  deviator  and  a  and  e  the  hydro- 
estatic  stress  and  strain  respectively.  The  kernels  ^(z)  and  K(z)  are  sums  of 
positive  decaying  exponential  functions,  i.e.. 


-X  z 

,  K(z)  =  I  K  e 
r 


(2.7a , b ) 


where  [J-r,  <*r,  Kr  and  X^  are  all  non-negative. 

In  the  generalized  endochronic  internal  variable  theory  the  evolution 
equations  are  expressed  in  terms  of  the  intrinsic  times  of  the  mechanisms  of 
internal  motion.  See  Ref  [22].  Specifically  to  each  q  the  theory  ascribes  an 
intrinsic  time  zr  such  that  the  equations  of  evolution  become 

dqr 

r  +  b  •  —rz —  =  0  (r  not  summed)  (2.8) 

.  r  r  az 
3q  -  r 


r  =  1  ,  2  . . .  n. 

In  the  endochronic  theory  of  plasticity  of  metal s  as  it  has  been  used  in  the 

past 

Z1  '  zZ  ‘  •••  ’  2„  '  2  ,2-?) 

and  the  elastic  bulk  modulus  K  in  constant,  i.e,  the  material  is  plastically 
incompressible.  Also, 


dz  =  ,  dc  =  i  deP i 


(2.10a ,b) 
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and 

P  d5 

dep  =  de  -  j—  (2.11 ) 

~  0 

where  Ug  is  the  elastic  shear  modulus.  The  function  f  is  positive  and  non 
decreasing.  Thus  in  the  case  of  plasticity,  and  in  the  context  of  the  above 
assumptions,  equations  (2.5)  and  (2.6)  have  the  form 


s  =  2u0(z)*de  (2.12) 

G  *  3Ke  (2.13) 

Substitution  of  equation  (2.11)  in  equation  (2.12)  gives  an  equivalent  consti- 
tuitive  equation  which  relates  s  directly  to  the  history  of  eP.  Thus 

s  =  2p (z)*deP  (2.14) 


See  Valanis,  Ref.  [23].  The  relation  between  p  and  is  given  in  terms  of  their 
Laplace  transforms  in  equation  (2.15). 

o'  (1  -  ^  )  =  v  (2.15) 

u0 

It  has  been  found  that  in  the  case  of  metals  at  room  temperature  p  and  f  are 
well  represented  by  the  relations 


1  -kz 

0  =  °o7  e 


f  =  1  -  ae 


■bz 


(2.16a  ,b) 


where  Pg,  a,  a,  b  are  positive  and  k  is  non-negative. 


In  the  Gibbs  formulation  the  thermodynamic  state  is  described  by  the  free 
energy  density  <t> ,  which  is  a  function  of  the  stress  tensor  a  the  temperature  T 
and  n  internal  variables  q  which,  again,  are  second  order  tensors.  The  func¬ 


tion  $  is  related  to  +  by  the  equation 
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*  =  U ')  -  0  .  .  E  .  . 

1J  U 


(2.16) 


The  counterparts  of  equations  (2.1)  and  (2.2)  are 


3  <J> 

e  *  '  To 


(2.17) 


(2.18) 


whereas  the  positive  rate  of  irreversible  entropy  gives  rise  to  the  inequality 


3<t>  *  r  *  r 

— -  •  q  >  0,  iq  n  *  0 

3qr  ~ 


(2.19) 


for  each  r,  where  -  is  the  internal  microforce  for  the  internal  mechanism  r 

-V 

r 

while  the  evolution  equations  for  q  are,  in  the  case  of  the  generalized 
endochronic  formulation 


3  4)  .  dqr 

+  b  •  -rr— 


r  dz 


(2.20) 


If  we  now  stipulate  that  4>  is  quadratic  in  its  variables  and  dz  =  dz  for 

r 

all  r,  and  br  are  constant  then  in  the  case  of  isotropic  materials  and  isothermal 

a 

conditions,  equations  (2.17)  and  (2.20)  combine  to  give  the  constitutive 


equations 


e  =  j  L(z)*ds 


(2.21) 


e  =  j  N(z)*do 


(2.22) 


When  plastic  incompressibility  applies  N  is  a  constant  and  equal  to  -  ( 

K 

equation  (2.13).  Also  u  and  L  are  related  by  equation  (2.23). 
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u(z)*dL  =  H ( z ) 

where  H(z)  is  the  unit  step  function. 
In  view  of  equation  (2.11) 


(2.23) 


pi5 
e  ~s'7u~ 
~  *  2  U0 


(2.24) 


ep  =  J(z)*ds 


(2.25) 


where 


J(z)  =  -  H(z)  +  Liz) 
u0 


(2.26) 


The  functions  J(z)  and  p(z)  are  also  related  as  shown  in  equation  (2.27) 


J(z)*dp  =  H(z) 


(2.27) 


Spectrum  of  Intrinsic  Times 


Of  the  n  internal  mechanisms  let  group  r  have  an  intrinsic  time  zr-  Then  to 

the  group  r  of  internal  variables  n^  in  number,  there  will  correspond  an  intrin- 

m 

sic  time  z  .  Evidently  J  n  *  n  where  m  is  the  number  of  groups. 
r  r=l  r 

A  straightforward  analysis  using  equations  (2.16)  (2.17)  and  (2.20)  when  4>  is 

(  p  ) 

a  quadratic  isotropic  function  of  e..  and  q..  and  b  are  constant  isotropic 

ij  „r 

tensors  leads  to  the  equations 


ep  =  l  J  (z  )*ds 
~  r=l  r  r 


(2.28) 


e=  0,  e  =  j  No 


(2.29a ,b) 


where  plastic  incompressibility  has  been  observed.  They  physical  1  nterpretation 
of  equation  (2.28)  is  that  each  of  the  m  groups  of  mechanisms,  say  r,  contributes 
to  the  total  strain  a  partial  strain  e£  such  that 


(2.30) 


where 

ej  -  OrUr)*ds 

***  ~ 


(2.31) 


The  intrinsic  time  zp  is  related  to  c  by  the  equation 


where  f  is  the  hardening  function  of  graph  r. 

The  need  for  this  more  general  approach  which  was  presented  in  a  previous 
reference  [22],  has  been  discussed  in  the  Introduction  and  haS  to  do  with  the 
fact  that  one  intrinsic  time  is  just  not  sufficient  to  describe  the  creep 
response  of  metals  to  piece-wise  constant  stress  histories.  The  physical  justi¬ 
fication  for  this  possibility  is  discussed  in  the  section  on  deformation  kine¬ 
tics,  but  simply  put,  it  means  that  the  hardening  function  f  is  not  the  same  for 
all  mechanisms  at  high  temperatures,  though  the  assumption  of  an  f  common  to  all 
qr  suffices  at  room  temperatures  as  demonstrated  in  our  work  on  endochronic 

mm 

plasticity.  The  appeal  to  deformation  kinetics  is  necessitated  by  the  desire  to 
determine  how  f  is  influenced  by  the  micro-mechanical  process  which  accompanies 


the  inelastic  deformation. 


3.  Deformation  Kinetics 


The  theory  of  deformation  kinetics  was  founded  by  Eyring  [24]  and  its  appli¬ 
cation  to  rate  processes  has  been  pursued  by  Eyring,  Krausz  [24]  and  their 
co-workers  with  a  great  deal  of  success. 

In  keeping  with  the  ideas  of  deformation  kinetics  we  attribute  macromotion 
to  additive  effects  of  micromotions  brought  about  by  local  distortion  of  atomic 
"energy  barriers".  At  this  point  it  is  essential  that  we  distinguish  between 
diffusion  of  particles  and  diffusion  of  dislocations  or  vacancies.  The  distinc¬ 
tion  may  be  stated  most  simply  in  terms  of  the  mean  free  path  I  of  a  particle. 

In  the  case  of  particle  diffusion  l  is  large  compared  to  its  counterpart  in  the 
case  of  dislocation  or  flaw  diffusion.  If  a  particle  travels  n  units  of  distance 
"a"  before  coming  to  rest,  then  l  =  na.  However,  if  a  flaw  travels  n  units  of 
distance  a,  the  mean  free  path  of  a  particle  is  still  a,  because  a  different 
particle  partakes  each  time  in  the  motion  of  the  flaw.  Thus  in  the  case  of  par¬ 
ticle  diffusion  n  >>  a.  Most  important,  however,  is  the  fact  that,  in  either 
case,  each  unit  of  motion  consists  of  an  atom  moving  across  an  energy  barrier. 

In  addition  we  would  expect  that  in  the  case  of  flaw  or  dislocation  diffusion 
the  energy  barriers  would  be  lower  than  those  of  particle  diffusion.  In  fact, 
the  activation  energy  of  self  diffusion  is  lower  at  low  homologous  temperatures 
(where  dislocation  motion  is  dominant)  than  at  higher  temperatures  where  par¬ 
ticle  diffusion  dominates  the  process. 

To  apply  the  ideas  of  deformation  kinetics  to  the  viscoplastic  deformation 
and  flow  of  solids  we  appeal  to  a  simple  atomic  model  whereby  prior  to  the 
application  of  stress  each  atom  of  the  solid  is  situated  at  the  bottom  of  a  sym¬ 
metric  potential  well.  A  typical  potential  well  with  the  accompanying  local 

potential  surface  is  shown  in  Figure  1  by  a  solid  line.  The  forward  and  back- 

|* 

ward  barriers  are  equal  and  both  have  a  height  . 


When  stress  is  applied, 
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the  atoms  will  be  displaced  from  their  initial  positions  and  the  local  potential 
surface  of  an  atom  will  distort.  The  distorted  potential  surface  is  also  shown. 
The  effect  of  the  distortion  is  to  reduce  the  forward  barrier  by  an  amount  w^ 

p 

and  increase  the  backward  barrier  by  an  amount  w^  .  It  will  be  shown  that  this 

type  of  barrier  distortion  will  give  rise  to  an  average  forward  motion  of  the 

r 

atoms  occupying  potential  wells  with  barriers  Eg. 


Boltzmann  Statistics.  To  determine  quantitatively  the  effect  of  barrier 
distortion  on  the  mean  atomic  motion  we  appeal  to  Boltzmann  statistics. 
Accordingly  the  probability  of  finding  an  atom  in  an  energy  state  e.  is  given  by 
equation  (3.1): 


-Be^ 

p^  =  ae  (3.1) 

where 

~®ei  1 

a  =  1/J  e  *  ®  =  TcT  (3.2a,b) 

i 


in  the  usual  notation. 

Let  Nr  be  the  number  of  atoms  occupying  potential  wells  with  initial  energy 

barriers  Eg  .  The  probability  that  an  atom  is  in  an  energy  state  greater  than 

r  .  r 
Eg  is  pQ  where 


Pn  =  [  a  exp  ( -6 e Q )  (3.3) 

r  r 

Ei>E0 

States  e^  such  that  eT  >  Eg  we  have  called  activated  states  [25]  differing  from 

Eyring.  Thus,  the  number  of  atoms  in  an  activated  state  is  N  p!^  .  As  the 

r  o 

barriers  are  symmetric  the  probability  that  an  atom  will  move  forward  is  equal 
to  the  probability  that  it  will  move  backwards  so  that  the  net  motion  (ave-age 
displacement)  of  the  atoms  Nr  is  zero. 
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When  the  potential  energy  surface  is  distorted  the  number  of  atoms  Ar  par¬ 
taking  in  a  forward  motion  is  now  changed  to 

A  =  N  I  a  exp  (-Be^)  (3.4) 

r  rr  r  r 

ei>e0'Wf 

while  the  number  B  of  atoms  partaking  in  as  backward  motion  is 
r 


Br  =  Nr  I  a  exp  ( -6 e T ) 


(3.5) 


r  r  r 

e  .  >e  -+w. 
i  0  b 


The  net  number  of  atoms  that  partake  in  a  forward  motion  is,  thus,  A  -  B  where 

r  r 


r  r  r 

e  .  <£  n+w. 

1  0  b 

B  =  N  l  a  exp  (-Be. ) 

r  rr  r  r  1 

e ■>£k"W< 
i  b  f 


(3.6) 


To  evaluate  the  sum  on  the  right  had  side  of  equation  (3.6)  we  shall  assume  that 
r  r 

and  wfc  are  both  small.  In  this  case  we  represent  the  distribution  of 
r  r 

energies  e i  in  the  vicinity  of  eQ  in  terms  of  the  local  tangent  to  the  distribu- 
.  r 

tion  at  by  writing 


*,r  •  *0  *  kr(1  ■  'o’ 


(3.7) 


where  ip  is  the  value  of  i  at  eT  =  and  kr  is  the  slope  of  the  distribution 

p 

which  is  a  function  (in  general)  of  .  See  Figure  2. 

Substitution  of  the  relation  (3.7)  in  the  sum  on  the  right  hand  side  of 


equation  (3.6)  leads  to  the  simple  expression 

r 


-Be' 
e  0 


A  -  B  =  2N  a 
r  r  r 


sinBrw 


1  -  e 


Bkr 


(3.8) 


where 
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W 


r 


r  r 
"f*b 


2 


.  r  V  wf 

Ae0  =  ~1~ 


(3.9,10) 


We  note  that  the  terms  e  0  and  1  -  e  are  related  to  the  initial  state,  i.e., 
the  barrier  height  and  the  energy  distribution,  while  the  terms  sin  Bwr  and 

exp  (-BAEq)  are  brought  about  by  barrier  distortion. 

The  mean  velocity  vp  relative  to  the  lattice  of  the  atoms  in  group  r  may  now 
r 

be  calculated  in  terms  of  the  barrier  distortion  parameters  Nr  and  Ae^  .  If 
x  is  the  mean  lattice  distance  and  t  is  the  average  time  taken  by  the  atoms  of 
group  r  to  traverse  that  distance  across  the  barrier  e^  then 

vP  =  (X/t  )  { A  -  B  )/N  (3.11) 

r  r  r  r  r 

We  now  define  an  internal  variable  q  by  the  relation 

r 


» 


q 


r 


Evidently  in  view  of  equation  (3.11) 


(3.12) 


q  =  (A  -  B  )/N  t 
r  r  r  r  r 

In  view  of  equations  (3.8)  and  (3.13) 
r 


-Be 


2a  e 


-BAe 


1 

q  =  - -  (— )e  sinh  B  w 

r  -Bk  r  r 


(3.13) 


(3.14) 


1  -  e 

In  so  far  as  steady  creep  is  concerned  the  assumption  is  usually  made  that 

of  all  the  operating  mechanisms  only  one  survives  in  the  steady  state,  i.e., 

r  =  1  and  w1  is  proportional  to  the  stress  (in  one-dimensional  stress  fields). 

However,  in  the  case  of  transient  creep  it  is  the  local  microforce,  i.e., 

-  j— -  on  the  group  r  that  will  determine  the  barrier  distortion.  Thus,  following 
qr  3(t 

Ref.  [25],  we  let  w  be  proportional  to  -  - —  according  to  equation  (3.15) 

r  d  Q 
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L 


W  =  -c  — 

Wr  r  3q 


(3.15) 


where  C  may  depend  on  temperature, 
r 

Thus  if  we  let 

r  -6kr  ScJ 

b00  *  (1  -  e  )  e  /2a 


(3.16) 


then  in  view  of  equations  (3.14  -  3.16) 

r 


.  r  *  1 

b00  «r‘r' 


-6Ag 


0  si nh  (B  Cp  |^~)  =  0 


(3.17) 


This  is  an  "internal  variable"  form  of  the  equation  for  the  mean  irreversible 
motion  of  the  group  r  of  paricles  facing  a  potential  barrier  . 


Discussion  of  equation  (3,17).  As  we  pointed  out  equation  (3.17) 

establishes  a  physical  meaning  for  the  internal  variables  in  that  qr  is  the  mean 

displacement  relative  to  the  lattice  of  a  group  r  of  particles  facing  a  poten- 

r 

tial  barrier  of  magnitude  eQ  .  The  above  equation  was  published  by  Valanis  and 

p  p 

Lalwani  in  Ref.  [25],  with  Agq  =  0.  The  appearance  of  the  term  Ae0  in  equation 
(3.1)  was  inferred  as  a  result  of  our  effort  to  describe  analytically  the  creep 
response  of  306  stainless  steel  to  piece-wise  constant  stress  histories.  This 
will  be  discussed  in  Section  4. 

The  time  to  traverse  the  barrier,  i.e.,  xr  is  also  of  central  importance  in 
equation  (3.17).  Eyring  used  simplifying  assumptions  to  arrive  at  the  conclusion 
that  is  proportional  to  the  square  root  of  the  ambient  temperature.  However, 
one  can  show  that  it  depends  at  least  in  part  on  the  barrier  shape  and  height 
(Ref.  [25]).  In  this  work  we  have  found  that  is  also  sensitive  to  the  plastic 
strain  rate.  This  is  to  be  expected  sind  tr  depends  on  the  barrier  confor¬ 
mation,  which  in  turn  depends  on  the  plastic  strain.  The  rate  of  plastic  strain 

affects  the  rate  of  barrier  distortion  which  must  affect  the  traversal  time  t  . 

r 
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Consider  now  two  processes  (a)  and  (b)  the  first  of  which  is  proceeding  at  a 
faster  plastic  strain  rate  than  the  second.  With  regard  to  the  forward  motion 
of  a  particle,  the  height  of  the  barrier  will  be  diminishing  faster  in  case  (a) 
than  in  (b),  so  that  the  forward  moving  particle  will  be  encountering  a  con¬ 
sistently  lower  barrier  in  case  (a)  than  in  (b).  It  follows  that  the  time  to 
cross  the  barrier  in  case  (a)  will  be  shorter  than  in  case  (b).  Thus 


r  <  tu  whenever  c  >  c. 

3D  3D 

The  above  inequality  will  be  satisfied  if 


g(i) 


(3.18) 


where  xQ  is  a  constant  and  g  is  a  monotonical ly  increasing  function  of  c.  In 
this  work  we  have  set 


9(c)  =  r,_m  (3.18a) 

where  m  is  a  material  constant. 

Deformation  kinetics  is  brought  into  accord  with  linear  irreversible  ther¬ 
modynamics  if  in  equation  (3.17)  the  argument  of  the  hyperbolic  sine  is  suf¬ 
ficiently  small  for  the  approximation 


sinh  (8  cr  }  ~  8  cr  ^ 


to  be  appropriate.  In  this  event  equation 
lution  equation,  i .e. , 


brq  +  =  0 

r  3qr 


(3.19) 

(3.17)  becomes  a  standard  linear  evo- 

(3.20) 


where 


BAc 

e 


r 

0 


(3.21) 


16 


In  view  of  equations  (3.18a),  (3.20)  and  (3.21)  one  finds  that  the  "endochroni c" 
form  of  equation  (3.20)  is 
r 


where 


r  dq  3 _  n 

b0  dT  *  3q;  -  ° 


sr  -  s  r  T° 

b0  "  b0O  BC 


(3.22) 


(3.23) 


and 


.  _  dz 


(3.24) 


BAe 


f  =  e 
r 


‘°  (3.25) 

Thus  in  deformation  kinetics  terms  the  rate  sensitivity  is  attributable  to  the 
time  to  cross  the  barrier  while  the  hardening  (softening)  is  related  to  the 
change  in  the  mean  height  of  the  barrier  as  a  result  of  the  stress  history. 

f* 

Thus  if  aeq  increases  the  material  hardens  while  if  it  decreases  the  material 
softens  in  accord  with  our  physical  intuition  regarding  such  processes. 


17 


4.  Analysis  of  Piece-wise  Constant  Stress  Histories  in  Pure  Shear. 


We  begin  with  the  integral 


eP  =  J(z)*ds  (4.1 ) 

where  eP  represents  a  shear  creep  strain  component,  s  is  the  corresponding  shear 
stress  component  and  J  the  appropriate  creep  function.  As  usual 


dz 


dc 

g(£)f 


(4.2) 


where  g  is  the  rate  sensitivity  function  and  f  the  hardening  function.  Also 


dc  =  k  |deP  |  (4.3) 

where  k  is  an  appropriate  scalar  constant.  Typically,  it  eP  denotes  a  creep 
shear  strain  component  and 

ds  =  >deP«  (4.4) 


then  k  =  /2  and  J  is  the  creep  function  in  pure  shear. 

For  our  purposes  it  is  more  convenient  to  write  equation  (4.1)  in  the  expli¬ 
cit  form 


t 

ep  =  /  J ( z ( t )  -  z ( t ' ) )  •rr-r  dt' 
0 

for  reasons  that  will  become  apparent. 


(4.5) 


4- .1 .  Monotonic  Creep  in  the  Presence  of  a  Constant  Stress  History 
In  this  specific  case 

s ( t )  =  sQ  H ( t)  (4.6) 

where  H(t)  is  the  unit  step  function  whose  "derivative"  is  the  Dirac  delta  func¬ 
tion.  In  this  instance  substitution  of  equation  (4.6)  in  equation  (4.5)  gives 
the  creep  response  in  the  simple  form 
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eP  .  s0  J(2)  U.7> 

where  $0  is  the  amplitude  of  the  step  function  of  applied  stress.  It  is 
apparent  from  equation  (4.7)  that  if  the  form  of  J(z)  is  known  then  knowledge  of 
z(t)  determines  creep  strain  in  terms  of  the  stress  amplitude  s„.  We  caution 
that  ep  is  not  necessarily  linear  in  s  since  z(t)  depends  on  sQ  as  we  shall 
demonstrate. 

To  this  end  differentiate  equation  (4.7)  with  respect  to  t  and  use  equations 
(4.2)  and  (4.3)  to  find  that  under  monotonic  conditions 

g(c)  =  ksQ  J'(z)f’ 1  (4.8) 

where  J'(z)  is  the  derivative  of  J  with  respect  to  z.  Thus 

5  =  9_1  (k  s0  J'(z)  f'1  }  (4.9) 

But  from  equations  (4.7)  and  (4.3) 


C  =  k  sQ  J 1 ( z )  z 

Thus  from  equations  (4.9)  and  (4.10) 


(4.10) 


.  g"1{ks0J,(z)f'1l 

2  =  k  SqJ'Iz) 

Equation  (4.11)  gives  z(t)  by  numerical  integration  if  f  is  known. 


(4.11) 


To  assign  analytical  forms  to  the  functions  J(z)  and  g(?)  we  appeal  to 
experiment  and  the  underlying  assumptions  of  endochronic  plasticity.  It  is  com¬ 
mon  experience  that  metals  become  more  strain  rate  sensitive  as  the  temperature 


rises.  However,  the  spirit  of  the  endochronic  theory  is  that  this  change  is 
brought  about  not  by  a  change  in  the  form  of  J(z)  but  by  virtue  of  g(;)  which 
is  evidently  dependent  on  temperature  even  though  this  dependence  is  suppressed 
in  equation  (4.2). 


At  room  temperature  where  rate  effects  are  not  significant  J(z)  is  repre¬ 
sented  very  closely  by  the  analytical  impression 

J(z)  =  J0  (4.12) 

where  a  is  the  vicinity  of  0.85  for  a  number  of  metals.  This  form  is  retained 
by  virtue  of  the  above  argument,  at  higher  temperatures. 

Experiments  also  indicate  that  under  monotonic  creep  conditions 

ep=F(s0)t6  (4.13) 

i.e.,  that  the  stress  and  time  dependence  of  creep  strain  are  factorable  and 
that  the  time  dependence  is  represented  very  closely  by  a  power  law.  If  during 
monotone  creep  f  is  a  constant  -  which  was  found  to  be  so  for  one  component  of 
the  creep  -  then  for  equation  (4.13)  to  hold  g(S)  must  also  be  a  power  function. 
Thus  we  have  set 

g(S)  =  im  (4.14) 

In  view  of  these  stipulations  equation  (4.11)  now  becomes 

2  l^m  (q-1 )(l-m) 

fm  2  -  (k  a  s0  JQ)  Z  115  (4.15) 

Special  solutions  to  equation  (4.15).  We  proceed  to  give  some  special  solutions 
to  equation  (4.15)  when  (a)  f  is  constant  and  (b)  when  f  is  a  power  function  of 
z.  When  f  is  constant  the  solution  is  given  by  equation  (4.16): 

il 

fm  z  =  A  (k  a  sQ  Jq)  t  8  (4.16) 

where 

A  =  (tt)8  ,  S’  =  m/1  +  am  -  a,  n’  =  S’  (4.17) 

p  rn 

In  the  case  where  f  =  fQ  z  ,  where  *  >  0,  the  solution  to  equation  (4.15)  is 
given  by  equation  (4.16)  as  before  except  that  now  the  constant  8’  is  given  by 
equation  (4.18): 


Discussion.  So  far  we  have  represented  the  creep  strain  by  a  single  integral. 

We  have  also  represented  the  creep  function  J(z),  the  strain  rate  sensitivity 
function  g(s)  and  the  hardening  function  f(;)  by  analytical  forms  of  the  power 
type.  By  analysis  we  then  arrived  at  equation  (4.19)  which  is  basically  of  the 
form 

eP  =  B  sj  tB  (4.21  ) 

where  B  is  a  constant,  whereby  the  monotonic  creep  strain  depends  multiplica- 
tively  on  the  stress  amplitude  to  the  power  n  and  the  time  to  the  power  B.  This 
form  has  appeared  frequently  in  the  literature  where  it  has  been  arrived  at  by 
analysis  of  the  data.  It  does  not  for  all  creep  data  and  certainly  not  over  the 
entire  range  of  stress. 

What  is  important,  however,  is  that  the  creep  strain  depends  on  time 
according  to  a  power  law  (equation  (4.21)  in  accordance  with  observation  as  per 
equation  (4.13)  while  the  dependence  on  stress  is  of  a  more  general  type.  One 
can  change  the  dependence  of  eP  on  s^  by  changing  the  analytical  form  of  J(z)  or 
g(S)  or  both  but  it  seems  that  this  would  vitiate  the  dependence  of  eP  on  a 
power  function  of  t.  Two  other  avenues  are,  however,  available.  One  is  to 
introduce  a  spectrum  of  intrinsic  times,  as  discussed  previously,  i.e.,  a  series 
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of  integrals  on  the  right  hand  side  of  equation  (4.1).  The  other  is  to  intro¬ 
duce  a  stress  dependence  in  the  hardening  function  f  .  This  has  been  found  to 
be  the  case  in  other  materials  such  as  polymers. 

Specifically,  if  one  sets 


f  ■  Wz 


(4.22) 


then  under  monotonic  creep  conditions  in  the  presence  of  constant  stress 
equation  (4.19)  will  have  the  form  given  by  equation  (4.13),  for  an  appropriate 
choice  of  the  function  f. 


This  approach  alone,  however,  has  been  found  inadequate  to  describe  creep 
under  cyclic  piece-wise  constant  stress  histories.  This  question  as  well  as  a 
constitutive  equation  involving  more  than  one  intrinsic  time  will  be  discussed 


in  the  next  section. 
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5 .  Specific  Constitutive  Equations  for  304  Stainless  Steel 

In  the  application  of  the  theory  to  304  stainless  steel  at  600°C  and  speci¬ 
fically  to  the  experimental  data  generated  by  Ohno  et  als.  (Ref.  [4]),  two 
terms  were  retained  on  the  right  hand  side  of  equation  (2.28),  i.e.. 


ep  =  J,(z,)*ds  +  JJzJ*ds 


(5.1) 


For  the  purposes  of  analysis  and  presentation  of  the  results  it  is  more  con¬ 
venient  to  write  equation  (5.1)  in  the  form 


V  _  P  ,  QP 
e  -  e,  +  e2 


(5.2) 


where 


;P  =  J  (z  )*ds 
r  r  r 


(5.3) 


r  =  1,  2.  In  this  case  two  hardening  functions  exist  in  the  sense  of  equation 


(5.4) 


,  dz  1 

dz_  =  7-  •  — 

r  fr  cm 


(5.4) 


where  m  is  a  material  constant  found  to  be  equal  to  0.12.  Also  two  creep  func¬ 
tions  J.|(z)  and  J ^(z)  are  needed  and  these  were  given  the  analytical  forms  shown 


in  equation  (5.5): 


J1  =  J1  2  1  ’  °2  =  °2  Z 


(5.5  a , b ) 


where  a  =  0.836,  a2  =  1 .  »  2.34  x  10"3  MPa,  =  1.58  x  10'3  MPa.  It  still 

remains  to  determine  the  form  of  the  hardening  functions.  In  order  to  match  the 
monotonic  data  was  represented  by  a  power  function  of  the  form 


f2  =  Z 


(5.6) 


where  =  0.196. 
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On  the  other  hand  f^  could  not  be  represented  as  a  state  function  of,  say,  s 
and  z,  or  any  other  variable  for  that  matter.  Rather,  the  experimental  cyclic 
results  of  Ref.  [4],  gave  strong  indication  that  f^  should  be  given  in  differen¬ 
tial  form  of  the  type 

d  log  f1  =  dF(  | s  | ,  z1 )  |2  (5.7) 

1 

Note  that  the  right  hand  side  of  equation  (5.7)  is  not  an  exact  differential  and 
hence  f  is  a  function  of  the  stress  history.  The  physical  implication  of 
equation  (5.7)  is  that  a  change  in  z.|  does  not  affect  ^  if  during  the  change 
the  absolute  value  of  the  stress  s  remains  constant.  A  mathematically  more 
explicit  form  for  f^  is 

d  log  f]  =  |^Tj-  (  |s  |,  2])  d  [ s  |  (5.8) 

The  logarithmic  form  is  not  fortuitous  but  is  a  consequence  of  the  physics  of 
deformation  kinetics  and  specifically  equation  (3.25)  in  view  of  which 


log  f]  =  B  A  c'Q  (5.9) 

d  log  f1  =  6  d  (A  z'Q)  (5.10) 


The  implication  is  that  in  mechanism  1  the  mean  barrier  height  will  change  when 
the  absolute  value  of  stress  changes  but  not  otherwise.  The  constitutive 
description  of  the  material  is  complete  once  the  function  F(|s|,  z^)  is  known. 
The  function  F  is  given  below  for  various  values- of  |s|  (in  MPa): 


F ( 1 34 . 2 ,  Z-| )  =  2.25  +  1.3(1  -  e'30zl  ) 
F ( 120,  z^j )  «  F ( 1 34 . 2 ,  z1 ) 

F(90,  z.| )  =  2.19  +  0.9(1  -  e‘5°21) 

F (60,  zn )  =  1.86  +  44(1  -  e"85zl) 


(5. 1  la ,b,c,d) 


In  addition  the  value  of  the  hardening  function  at  zero  stress  and  zero  value  of 
z  is  set  at  0.1054.  This  value  together  with  the  relation. 

log  f  (0,  |s  |n )  -  log  f{0,js|2)  =  FfOjs^  -  F(0,|s|2)  (5.12) 

where  s^  and  s2  are  any  two  stress  levels  determines  f  for  various  values  of  the 
initial  stress  applied  at  the  onset  of  a  creep  experiment. 
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Monotonic  Creep  Experiments 

The  strain  response  to  a  constant  stress  history  is  obtained  by  application 
of  equation  (4.19)  in  conjunction  with  equation  (5.2).  Specifically 


P  ?  ,  -{Pr 

e  =r=l‘  0r 


®  f* 

"“r1  (SoV  ‘ 


(6.1) 


in  the  presence  of  the  constraint 


6 i  =62=6,  n1  =n2=n,  m]  =m2=m 
so  that  equation  (6.1)  becomes 


=P  ■  <V"  ‘‘  I,  *  r«V  "  “V"''  V 

r*  I 


(6.2) 


(6.3) 


In  the  case  of  the  linear  model  (2),  fQ2  is  a  constant.  However,  in  the 
non-linear  model  (1),  fQ1  is  a  function  of  |s|.  This  dependence  is  determined 
by  adjusting  fg,  for  various  values  of  |s|,  so  as  to  obtain  optimal  agreement 
between  theory  and  experiment  in  the  case  of  monotonic  creep.  The  function 
f  (  |S  |)  is  shown  in  Figure  3.  With  all  the  other  constants  known,  the 
descriptive  capability  of  the  theory  is  shown  in  Figure  4. 


Cyclic  Creep  Experiments 

Following  Section  2  let  e^  be  the  r'th  partial  shear  strain  such  that 


e„  =  J  (z  )*ds  (6.4) 

r  r  r 

and 

ep-?ef  (6.5) 

r=l 

where  in  our  case  m  =  2.  For  our  purposes  it  is  more  convenient  to  write  equation 
(6.4)  in  the  form 
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eP  =  /  JrCzr( t)  -  z(tr)]  £fv  d t 1 


(6.6) 


in  the  specific  case  of  piece-wise  constant  stress  histories  of  the  type  con¬ 
sidered  by  Ohno  et  als. 


#  =  sj5(t)  -  a  6(t-t.)  +  a  S(t-tJ  ....  } 
at  u  I  ^ 


(6.7) 


where  6 ( t)  is  the  Dirac  delta-function,  sn  is  the  initial  stress  amplitude  and  a 
is  a  constant.  In  this  set  of  experiments  two  parameters  sQ  and  "a"  define  the 
history  of  stress  -  in  addition  to  the  reversal  times  t^ ,  t^  ...  t^. 

Substitution  of  equation  (6.7)  in  equation  (6.6)  and  integration  gives  the 
explicit  result 


>p  =  S0  Urur)  +  a  l  <-l)r  Jr(zr-2r„)| 
n=  1 


(6.8) 


'  s0  ,Ur; 


(6.9) 


where  y(zr)  represents  the  bracket  on  the  right  hand  side  of  equation  (6.8).  We 


differentiate  equation  (6.9)  to  obtain 


K  -  s0  r 


(6.10) 


where  y  =  .  Use  of  equation  (5.4)  then  gives  the  result 

r  dzr 


eP  =  s  y 
r  S0^frZ"« 


(6.11) 


Now  we  take  absolute  values  of  both  sides  of  equation  (6.11)  and  use  equation 
(4.3)  to  obtain 
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•  m 

fr'  =sok|yl 
But  in  view  of  equation  (5.4 

1 

1  -m 


(6.12) 


(6.13) 

Equations  (6.12)  and  (6.13)  combine  to  give  the  following  differential  relation 


4  =  (f  ZJ 
r  r 


between  dz  and  dt 
r  1 


dt  = 


f  m  . 

r  dzr 
7,  -  ,1 
'k$0  yr!-i" 


(6.14] 


Integration  of  equation  (6.14)  gives  the  relation  between  zp  and  t. 

Substitution  of  z  (t)  in  equation  (6.8)  then  gives  the  desired  relation  ep(t) 
r  r 

and,  therefore,  eP(t)  upon  use  of  equation  (6.5). 


In  our  particular  case 


Jr  ~  °0r  Zr 


and 


al  f'  r  -r 

a  j  z  +  al  (-1 )  Jn  (a  -  a  ) 
r  Or  r  Or  r  rn 


(6.15) 


(6.16) 


In  Figures  5,  6  and  7  we  show  the  experimental  values  of  eP(t)  obtained  in  Ref. 
[4]  for  the  stress  histories  shown.  Also  shown  are  the  analytical  functions 
eP(t)  obtained  (a)  by  the  use  of  the  present  theory  and  (b)  as  reported  by  Ohno 
et  als.  in  Ref.  [4]  using  their  own  theory. 

The  most  significant  difference  in  the  predictive  capability  of  the  two 
theories  lies  in  their  depiction  of  the  creep  recovery  slope  at  points  of  stress 
reversal.  The  endochronic  theory  predicts  an  infinite  slope,  in  agreement  with 
experiment,  while  the  theory  by  Ohno  et  als.  depicts  a  finite  much  shallower 
slope.  Overall  the  predictive  capability  of  the  endochronic  theory  is  very  good 
for  the  type  of  stress  histories  discussed  here. 


Creeo  in  Tension-Tors io 


In  this  section,  we  examine  this  problem  in  the  complex  case 
where  both  the  torsion  and  tension  histories  are  piece  wise 
constant  junctions  of  time.  In  addition  in  the  complexity  of  the 
stress  histories,  we  have  the  added  coupling  effect  which  is 
observed  experimentally.  In  other  words,  the  cyclic  creep  in 
torsion  is  affected  by  the  presence  of  tension.  It  will  show  that 
this  effect  is  accounted  for  satisfactorily  by  the  definition  of 
intrinsic  time.  In  this  section  we  shall  address  the  experiments 
of  ohms  [4]  where  cyclic  creep  in  shear  is  carried  out  at  constant 
tensile  stress. 


The  basic  equations  are; 


Bb  -  i  «  r 
r=l 


(7.1) 


f  =  £(Zr)  * 


(7.2) 


fr  * 


(7.3) 


1 1  or  i 


(7.4) 


Jr  Z 


(7.5) 


(7.6) 


Thus  the  actual  solid  is  represented  by  two  inelastic  model  solids 
in  series.  The  pertinent  material  functions  and  constants  are 
the  following: 


al  * 

0.836, 

OrH 

hj 

*  2.34X10"3 

MPa 

a  2  " 

1.0 

J2 

*  1.58X10-3 

MPa 

®1  * 

m2  -  0 

.12 

*1  * 

o 

=*  0 . 
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In  this  section  f^  is  a  constant  for  simplicity  while  in  section  6 
it  was  an  elaborate  function  designed  to  match  the  experimental 
data  as  closely  as  possible.  Thus  in  this  section,  the 
effectiveness  of  the  theory  is  demonstrated  without  the  need  for 
the  complexity  of  representation  of  f^  adopted  in  section  6. 

The  stress  histories  involved  in  the  experiments  of  Ref.  4  are 
of  the  type 

aij(k)  -  4j  Ha(t)  (7.7) 

where 

Ha(t)  -  H (t )  -  a  £>1  H(t-rt0)  (7.8) 

and  H ( t)  denotes  the  Heaviside  step  function. 

In  the  following  denote  the  axial  stress  by  aa,  the  axial 
creep  strain  by  eP,  creep  strain  by  eP.  In  all  the  experiments 
considered  below 

al  •  al  H(t) 


(7.9) 


while 


(7  .  10) 


The  creep  strains  of  eP  and  eP  were  then  determined  numerically 
following  the  method  given  in  detail  in  section  6.  The  results 
are  shown  in  Figures  8,  9  and  10  for  corresponding  values  of  a=l, 
1.5  and  2.  Comparison  with  data  shows  good  qualitative  agreement. 
However  a  refinement  of  the  representation  of  the  hardening 
properties  of  the  material  seems  pertinent. 

The  work  of  this  section  was  done  co-jointly  by  K.C.  Valanis 
and  S.D.  Lee.  The  latter  is  a  graduate  student  in  Civil 
Engineering  at  the  University  of  Cincinnati. 


In  this  section  we  extend  the  theory  to  the  domain  of  large  deformation.  This 
is  am  important  development  because  high  tenperature  processes,  during  which 
materials  are  strain  rate  sensitive,  involve  large  deformations,  air  target 
here  is  metals. 

A  great  deal  of  effort  has  been  expended  in  recent  years  to  develop 
constitutive  eqatians  pertaining  to  large  deformation  of  metals.  The  problems 
that  one  encounters  in  the  formulation  of  small  deformation  theories  are 
magnified  when  large  deformations  are  involved.  When  the  problem  is  approached 
from  a  yield  surface  point  of  view,  one  has  to  be  concerned  with  the  evolution 
of  its  geometry  in  stress  space  as  the  material  deforms  and  its  translation  in 
stress  space  (  i.e. ,  a  constitutive  description  of  the  back  stress  -  which  has 
been  the  central  problem  in  classical  plasticity) . 

In  addition,  and  irrespective  of  one's  approach,  the  question  of 
appropriate  separation  of  the  increment  of  plastic  strain  (defined  in  large 
deformation  terms)  into  elastic  and  plastic  parts  must  be  addressed  and 
resolved.  Hew  this  is  to  be  done  is  a  matter  of  differing  opinion.  Of  course 
over  and  above  these  considerations,  the  principle  of  isotropy  of  space  must 
must  apply.  This  last  requirement  creates  other  difficulties  associated  with 
"objective  rates"  when  incremental  theories  are  considered.  In  strictly 
mathematical  theories  the  choice  of  a  physically  "correct"  objective  stress 
rate  (say)  is  not  obvious  a  priori. 

In  the  present  paper  we  side-step  a  number  of  the  above  problems  by 
utilizing  the  concepts  of  endochronic  plasticity  in  can junction  with  the  theory 
of  internal  variables  to  arrive  at  a  constitutive  equation  of  the  hereditary 
type  whereby  the  stress  tensor  (in  terms  of  its  covariant  components  in  the 
material  frame)  is  a  quadratic  functional  of  the  history  of  the  Right  Cauchy- 
Green  tensor.  Plastic  incompressibility  is  assumed  and  the  elastic  deformation 
is  neglected  relative  to  the  large  plastic  deformations  considered  in  the 
paper. 


The  theory  is  applicable  to  both  strain  rate  dependent  and  strain  rate 
independent  processes.  In  the  applications  a  constant  strain  rate  has  been 
assumed.  This  is  equivalent  to  a  constant  plastic  strain  rate,  in  the  light  of 
the  assumption  that  the  elastic  ceuponerrt  of  strain  is  negligible. 

A  number  of  solutions  of  practical  interest  have  been  obtained  in 
closed  form  or  by  finite  element  techniques  in  strictly  Lagrangian  terms  thus 
obviating  the  numerical  difficulties  associated  with  Euler ian  formulations. 
Also  difficulties  associated  with  objective  rates  are  not  encountered  because 
of  the  integral  form  of  the  constitutive  equation.  Solutions  associated  with 
bending,  torsion  and  inflation  sue  obtained  in  closed  form.  Finite  element 
solutions  pertaining  to  forging  of  blocks  have  also  been  obtained  and  will  be 
discussed. 

The  constitutive  equation  was  first  tested  by  application  to  the 
simple  homogeneous  extension  of  a  bar  in  the  presence  of  a  memory  kernel  that 
gives  rise  to  asynptotically  constant  Cauchy  stress.  Monotonic  behavior  was 
observed  with  no  instabilities.  The  Piola  stress  first  increased  and  then 
decreased  with  stretch  as  observed  in  experiment.  The  problem  of  plastic  flow 
in  a  tube  (extrusion  problem)  was  then  solved  in  closed  form.  Substantially 
flat  displacement  profiles  were  obtained,  in  agreement  with  observed  behavior. 

The  problem  of  the  large  inflation  of  a  thick  sphere  was  then 
solved  again  in  closed  form  revealing  a  geometric  instability  at  a  critical 
value  of  the  internal  pressure,  as  is  commonly  observed.  Finally  the  finite 
bending  of  a  beam  was  solved  revealing  the  observed  shift  of  the  neutral  axis 
toward  the  ccnpressive  side  and  strongly  ncn-linear  stress  distribution  within 
the  beam. 

The  "upsetting"  of  a  block,  i.e.,  forging  by  means  of  a  vertically 
applied  displacement,  was  solved  by  finite  element  methods.  The  initial 
barrelling  eventually  gave  way  to  a  bone-shaped  configuration  and  the  vertical 
stress  at  the  outside  boundary  went  from  ccnpressive  to  tensile  as  expected. 
The  ease  of  the  computation  is  enphasized. 

The  body  of  this  work,  which  served  as  the  Fh.D.  thesis  of  Dr.  J. 


Wang,  is  not  given  in  this  section,  but  is  appended  as  Appendix  I. 


References 


1.  Valanis,  K.C.,  "Endochronic  Theory  with  Proper  Hysteresis  Closure 
Properties,"  Systems,  Science  and  Software,  LaJolla,  CA,  EPRI  Report 
NP-1388,  1980. 

2.  Valanis,  K.C.,  and  Lee,  C.F.,  "Some  Recent  Developments  of  the 
Endochronic  Theory  with  Applications,"  Nuclear  Engineering  and  Design, 
69,  pp.  327-344,  (1982). 

3.  Valanis,  K.C.,  and  Lee,  C.F.,  "Endochronic  Theory  of  Cyclic  Plasticity 
with  Applications,"  J.  of  Appl.  Mech.  Trans.  ASME  (to  appear). 

4.  Ohno,  J.,  Murakami,  S.  and  Veno,  T.,  "A  Constitutive  Model  of  Creep 
Describing  Creep  Recovery  and  Material  Softening  Caused  by  Stress 
Reversals,"  Journal  cf  Eng.  Mat.  and  Technol.,  107,  1,  (1985). 

5.  Chaboche,  J.L.,  Dang  Van,  K.,  and  Cordier,  K.,  "Model i zation  of  the 
Strain  Memory  Effect  on  the  Cyclic  Hardening  of  316  Stainless  Steel," 
Trans.  5th  SMIRT,  Lll/3,  Aug.  1979,  Berlin,  Germany. 

6.  Krieg,  R.D.,  Swearengen,  J.C.,  and  Rohde,  R.W.,  "A  Physically-Based 
Internal  Variable  Model  for  Rate — Dependent  Plasticity,"  in  Inelastic 
Behavior  of  Pressure  Vessels  and  Piping  Components,  ed.  T.Y.  Chang  and 
E.  Krempl ,  ASME  PVP-028,  pp.  15-28,  (1978). 

7.  Malvern,  L.E.,  "Tne  Propagation  of  Longitudinal  Waves  of  Plastic 
Deformation  in  a  Bar  of  Material  Exhibiting  a  Strain-Rate  Effect,"  J. 
Appl.  Mech.,  Jjl,  pp.  203-208,  (1951). 

8.  Haisler,  W.,  "Application  of  an  Uncoupled  Elastic — Plastic  Creek 
Constitutive  Model  to  Metals  at  High  Temperatures,"  in  NASA  Symposium 
on  Nonlinear  Constitutive  Relations  for  High  Temperature  Applications, 
May  19-20,  (1982).  The  University  of  Akron,  Akron,  Ohio,  pp.  185-189. 

9.  Bradley,  W.L.,  and  Yuen,  S.,  "A  New  Coupled  Viscoplastic  Constitutive 
Model,"  in  NASA  Symposium  on  Nonlinear  Constitutive  Relations  for  High 
Temperature  Applications,  May  19-20,  (1982).  The  University  of  Akron, 
Akron,  Ohio,  pp.  217-233. 

10.  Cescotto,  S.,  and  Leckie,  F.A.,  "Determination  of  Unified  Constitutive 
Equations  for  Metals  at  High  Temperature,"  Proceedings  Int.  Conf.  on 
Constitutive  Laws  for  Engineering  Materials,  ed.  C.S.  Desai  and  R.H. 
Gallagher,  Jan.  10-14,  (1983),  Tucson,  Arizona,  pp.  105-111. 


11.  Krempl,  E.,  "The  Role  of  Servocontrolled  Testing  in  the  Development  of 
the  Theory  of  Viscoplasticity  Based  on  Total  Strain  and  Overstress," 
ASTM  STP  765,  pp.  5-28,  (1982). 

12.  Walker,  K.P.,  "Research  and  Development  Program  for  Nonlinear 

Structural  Modeling  with  Advanced  Time— Temperature  Dependent 

Constitutive  Relationships,"  NASA  CR-165533,  Nov.  1981. 

13.  Miller,  A.K.,  "An  Inelastic  Constitutive  Model  for  Monotonic  Cyclic, 
and  Creep  Deformation:  Part  I— Equations  Development  and  Analytical 
Procedures,  Part  II — Application  to  Type  304  Stainless  Steel,"  J.  of 
Eng.  Mat.  Techn.  Trans.  ASME,  98,  pp.  97-133,  (1976). 

14.  Hart,  E.W.,  "Constitutive  Relations  for  the  Nonelastic  Deformation  of 
Metals,"  J.  Engn.  Mat.  Techn.  Trans.  ASME,  98,  pp.  193-202,  (1976). 

15.  Ohashi,  Y.,  Ohno,  N.  and  Kawai ,  K.,  "Evaluation  of  Creep  Constitutive 
Equations  for  Type  304  Stainless  Steel  Under  Repeated  Multiaxial 
Loading,"  J.  Eng.  Mat.  and  Techno!.,  104,  159,  (1982). 

16.  Murakami,  S.  and  Ohno,  N.,  "A  Constitutive  Equation  of  Creep  Based  on 
the  Concept  of  a  Creep  Hardening  Surface,"  Int.  J.  Solids  and  Struct., 
J8,  597,  (1982). 

17.  Krausz,  A.S.  and  Faucher,  B.,  "A  Kinetics  Approach  to  the  Derivation 
and  Measurement  of  the  Constitutive  Equations  of  Time  Dependent 
Deformation,"  STP  765,  ASTM  (1982). 

18.  Wu,  H.C.  and  Yip,  M.C.,  "Strain  Rate  and  Strain  Rate  History  Effects 
on  the  Dynamic  Behavior  of  Metallic  Materials,"  Int.  J.  Solids  and 
Structures,  ]_6,  pp.  515-536,  (1980). 

19.  Wu,  H.C.  and  Yip,  M.C.,  "Endochronic  Description  of  Cyclic  Hardening 
Behavior  for  Metallic  Materials,"  J.  Appl.  Mech.  Trans.  ASME,  pp. 
212-217,  (1981). 

20.  Lin,  H.C.  and  Wu,  H.C.,  "On  the  Rate-Dependent  Endochronic  Theory  of 
Viscoplasticity  and  its  Application  to  Plastic-Wave  Propagation,"  Int. 
J.  Solids  and  Structures,  _1_9,  pp.  587-599,  (  1983). 

22.  Valanis,  K.C.,  "Proper  Tensorial  Formulation  of  the  Internal  Variable 
Theory,"  Archives  of  Mechanics,  29,  173,  (1977). 

23.  Valanis,  K.C.,  "Fundamental  Consequences  of  a  New  Intrinsic  Time 
Measure:  Plasticity  as  a  Limit  of  the  Endochronic  Theory,"  Arch,  of 
Mech.,  32,  pp.  171-191,  (1980). 

24.  Kraus,  A.S.  and  Eyring,  H.,  "Deformation  Kinetics,"  Wiley,  N.Y. 
(1975). 

25.  Valanis,  K.C.  and  Lalwani,  S.,  "Thermodynamics  of  Internal  Variables 
in  the  Content  of  the  Absolute  Reaction  Rate  Theory,"  J.  Chem.  Phys., 
67,  3980,  (1977). 


Fig.  2  Energy  Distribution  in  Vicinity  of  e0 
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Fig.  5  Creep  Strain  Response  to  Piece-wise  Constant 
Shear  Stress  History  (/5S  =  120  Mpa,  a=  1.0) 
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ABSTRACT 


The  object  of  this  research  is  an  extensive  study  of  an  endochronic 
constitutive  relation-  for  both  incompressible  and  compressible  materials 
undergoing  finite  plastic  deformation.  Attention,  in  this  dissertation,  has 
been  paid  to  the  foundation,  application,  and  computational  capability  of 
the  theory. 

In  the  first  part  of  this  dissertation  the  endochronic  constitutive 
relation  of  plasticity  for  incompressible  finite  deformation,  as  proposed  by 
Valanis  in  1978,  is  reviewed.  This  relation  is  formulated  in  the  material 
frame  of  reference  which  is  convenient  for  problems  Involving  a  definite 
initial  configuration  such  as  a  solid  as  opposed  to  a  fluid.  A  set  of 
academic  and  practical  problems  of  interest,  consisting  of  finite  plastic 
uniform  extension,  finite  plastic  shear  flow  in  a  pipe,  finite  plastic 
torsion,  finite  plastic  bending,  and  finite  plastic  spherical  expansion,  are 
solved  analytically.  Closed  form  solutions  are  obtained  for  all  problems  by 
the  use  of  a  semi-inverse  method. 

In  the  second  part  of  this  dissertation  the  theory  is  extended  to 
account  for  plastic  compressibility,  incompressibility  is  pivotal  in  the 
development  of  classical  theory  of  plasticity.  However,  it  is  only  an 
assumption  and  a  simplification  of  the  actual  situation.  In  the  endochronic 
theory,  plastic  compressibility  can  be  accommodated.  This  is  done  by 
modifying  the  free  energy  density  function  of  the  plastic  deformation 
process  by  adding  a  term  ip 0  ( I , )  which  reflects  the  compressible  plastic 
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deformation  and  then  developing  the  constitutive  relation  for  a  special 
functional  form  of  ^0(I,). 


To  solve  more  general  engineering  problems,  a  numerical  method  is 
developed  using  the  powerful  finite  element  technique  for  boundary  value 
problems  for  both  compressible  and  incompressible  materials.  In  the 
compressible  case,  a  special  form  of  ij;0  is  used  to  demonstrate  the 

application  of  the  theory.  The  finite  element  formulation  is  referred  to  the 
material  system,  by  using  a  Lagrangian  formulation. 

A  computer  code  is  then  established  to  solve  a  plane  strain  boundary 
value  problem  in  the  presence  of  compressible  plastic  deformation  by  the  use 
of  linear  triangular  elements.  A  specific  problem  associated  with  a  metal 
forging  process,  that  of  "upsetting"  a  block  is  analyzed  numerically  using 
the  code.  All  relevant  parameters  of  the  problem  are  investigated.  The 
results  obtained  give  a  very  reasonable  description  of  the  forging  process. 

The  study  of  this  dissertation  shows  that  the  endochronic  theory,  which 
is  based  on  a  sound  thermodynamic  foundation,  also  has  a  powerful 
computational  capability  for  solving  practical  engineering  problems  that 
involve  large  plastic  deformation. 
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CHAPTER  1  -  INTRODUCTION 


Finite  plasticity  as  a  subject  dealing  with  a  time-independent,  rate- 
independent,  large  permanent  strain  during  a  deformation  process,  has  very 
significant  application  in  engineering  problems.  In  ductile  metals,  under 
favorable  conditions,  plastic  deformation  can  continue  to  a  very  large 
extent  without  failure.  For  instance,  except  for  castings,  which  are  formed 
from  the  liquid  state,  all  metal  products  are  subjected  to  at  least  one 
metal  forming  process  during  their  manufacture.  In  metal  forming  processes 
such  as  forging,  drawing,  extrusion,  rolling,  stamping  and  cutting,  etc., 
the  products  suffer  a  considerable  3hape  change.  The  deformation  is 
substantially  permanent  and  involves  predominantly  large  plastic  strains. 
Other  processes  include  the  plastic  bending  of  beams  and  plates,  the 
overstraining  of  spheres  and  cylinders  which  are  widely  used  as  pressure 
vessels  in  the  chemical  industry  for  example.  A  thorough  understanding  of 
the  mechanics  involved  in  large  plastic  deformation  is  very  important  to 
engineering  application  and  design.  Moreover,  the  advancement  of  many 
branches  of  solid  mechanics  such  as  fracture  and  fatigue  as  well  as  soil 
mechanics,  rock  mechanics,  geophysics,  and  geology,  etc.,  is  closely  related 
to  the  development  of  a  sound  plasticity  theory.  Recently  with  the 
availability  of  high-speed  computing  facilities  and  the  development  of  the 
sophisticated  material  testing  machines,  many  researchers  have  been 
motivated  to  develop  more  advanced  theories  of  plasticity.  Hence  many  of  the 
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simplifying  assumptions  in  plasticity  and  many  empirical  formulations  are  no 
longer  necessary. 

A  difficult  element  of  classical  plasticity  is  the  concept  of  the  yield 
surface,  which  leads  to  experimental  and  numerical  difficulties  in  attempts 
to  describe  and  analyze  the  two-  or  three-  dimensional  response  of  a 
material.  Valanls  [1,2]  circumvented  this  difficulty  by  proposing,  in  1971, 
the  endochronlc  theory,  which  does  not  require  the  concept  of  yield  for  the 
description  of  plastic  behavior  of  materials.  During  about  15  years' 
development,  a  number  of  publications  documented  the  potential  of  the  theory 
to  describe  the  mechanical  response  of  materials  under  conditions  of  small 
plastic  deformation  and  proved  that  the  endochronlc  theory  of  plasticity  was 
able  to  predict  not  only  the  salient  features  of  the  plastic  behavior  of 
materials,  but  also  a  number  of  observed  features  of  plasticity  that  lay 
beyond  the  scope  of  the  existing  plasticity  theories.  It  is  natural  then  to 
extend  the  endochronlc  theory  to  the  other  subjects  in  plasticity  such  as 
problems  of  large  deformation,  rate-dependence,  and  thermomechanical 
coupling,  etc.  In  this  dissertation,  we  study  extensively  how  the 
endochronlc  theory  earn  be  applied  to  problems  in  the  domain  of  large  plastic 
deformation  and  how  it  can  be  used  to  solve  relevant  engineering  problems. 
Before  proceeding  with  the  development  of  the  theory,  we  will  review  briefly 
the  history  of  the  classical  and  endochronlc  theories  of  plasticity. 

Plasticity  as  a  science  is  generally  regarded  to  have  begun  in  186^ 
when  Tresca  C3]  published  a  preliminary  account  of  his  experimental  results 
on  punching  and  extrusion  and  formulated  a  yield  criterion  which  states 
that  a  metal  yields  plastically  when  the  maximum  shear  stress  attains  a 
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critical  value.  Early  contributions  to  the  theory  of  plasticity  were  also 
due  to  Saint-Venant  [4]  and  Levy  [5],  who  in  1870  applied  the  Tresca  yield 
criterion  to  establish  relations  between  stress  and  plastic  strain-rate  for 
two-dimensional  and  three-dimensional  plastic  deformation. 

In  the  following  60  years  development  of  the  theory  of  plasticity  was 
slow.  After  1921,  there  were  important  contributions  from  Von  Mises  [6], 
Hencky  [7],  and  Prandtl  [8].  Von  Mises  suggested  a  yield  criterion  on  the 
basis  of  purely  mathematical  considerations  and  Hencky  later  interpreted 
that  this  yield  criterion  implies  that  yielding  occurs  when  the  elastic 
shear-strain  energy  reached  a  critical  value.  Prandtl  showed  the  plane 
plastic  strain  problem  was  hyperbolic  and  calculated  the  loads  needed  to 
indent  a  planar  surface.  Hencky  continued  Prandtl's  work  and  discovered 
simple  geometrical  properties  of  the  field  of  slip-lines  in  a  state  of  plane 
plastic  strain.  Lode  [9]  and  Taylor  and  Quinney  [10]  carried  out  experiments 
on  various  metals  under  combined  tension  and  internal  pressure.  The 
effective  application  of  plasticity  theory  to  technological  processes  began 
in  1  925  when  Von  Karman  [11]  analysed  the  stress  distributior»  during  the 
rolling  of  metal  strip  by  an  elementary  method.  In  the  following  year  Siebel 
[12]  and  Sachs  [13]  presented  similar  theories  for  wire  drawing. 

About  1950,  the  classical  mathematical  theory  of  plasticity  entered  a 
fully  developed  period.  D.C.  Drucker  [14]  generalized  the  meaning  of  work 
hardening  and  related  it  to  the  stability  of  plastic  deformation  by 
postulating  that  (1)  the  work  done  by  an  external  agency  during  the 
application  of  additional  stresses  is  positive  for  a  working-hardening 
material,  and  (2)  the  net  work  done  by  an  external  agency  during  a  cycle  of 
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addition  and  removal  of  stresses  in  a  material  undergoing  plastic 
deformation  is  positive.  This  led  to  the  normality  of  plastic  strain  rata 
vector  with  respect  to  the  yield  surface  and  convexity  of  the  yield  surface. 
The  Von  Mises  flow  rule  was  then  derived,  which  provided  the  relationship 
between  the  loading  (yielding)  function  and  plastic  flow.  In  1961,  Ilyushin 
[15]  stated  his  postulate  of  plasticity,  which  he  claimed  to  be  a 
generalization  of  Drucker's  postulate.  A  third  aspect  of  the  classical 
plastic  theory,  besides  the  initial  yield  criteria  and  an  associated  flow 
rule,  Includes  the  isotropic  hardening  rule  by  Hill  [16]  and  Hodger  [17], 
kinematic  hardening  rule  by  Prager  [18],  modified  kinematic  hardening  rule 
by  Ziegler  [19],  and  Mroz’s  rule  of  hardening  moduli  [20]. 

The  theories  to  describe  the  relation  between  stresses  and  strains  are 
of  two  general  classes,  which  are  total  strain  (deformation)  theory  and 
incremental  strain  (flow)  theory.  The  total  deformation  theory  is  not 
physically  sound  because  it  cannot  account  for  history  effects  in  the 
mechanical  response  of  dissipative  materials.  Thus,  preference  is  given  to 
the  flow  theory,  which  is  still  useful  in  some  problems  [21,22].  Based  on 
the  classical  theory,  a  number  of  researchers  proposed  different 
constitutive  relations  for  large  plastic  deformation  problems.  They  include 
Hill  [23.24],  Rice  [25,26],  Mandel  [27,28],  Lee  [29,30],  and  Green  and 
Naghdi  [31].  However,  the  existing  results  in  the  literature  showed 
anomalous  solutions  to  analysis  of  necking  and  localization  in  metal  and 
peculiar  results  for  predic  ion  of  oscillatory  shear  stress  due  to  a 
monotonlcally  increasing  simple  shear  strain  [32-35]  when  the  kinematic 
hardending  rules  are  used. 


In  the  late  1960a,  the  formulation  of  constitutive  theories  of 
viscoelastic  materials  from  concepts  of  irreversible  thermodynamics  and 
Internal  state  variables  reached  an  advanced  level  of  development  [36].  On 
the  basis  of  this  success,  the  theory  of  plasticity  was  re-examined  along 
the  lines  of  irreversible  thermodynamics,  since  by  nature  plastic 
deformation  should  be  considered  as  an  Irreversible  thermodynamic  process. 
In  1971,  Valanis  [1,2]  proposed  a  new  approach,  called  endochronic  theory, 
for  describing  the  behavior  of  viscoplastic  material.  The  theory  was 
developed  on  the  concept  of  stress  and  free  energy  being  a  functional  of  the 
entire  histories  of  deformation  and  temperature  based  on  thermodynamics  and 
Internal  variables.  Two  new  concepts  were  drawn  in  the  development  of  the 
theory,  which  are  (1)  the  concept  of  intrinsic  time,  defined  as  the  norm  of 
the  Increment  of  the  total  strain  tensor.  It  is  a  scale  with  respect  to 
which  the  memory  of  a  material  of  its  past  deformation  history  can  be 
measured  and  (2)  the  concept  of  the  description  of  plasticity  without  a 
yield  surface.  The  application  of  the  theory  was  given  by  Valanis  [2,37].  Wu 
»nd  Lin  [38]  used  the  theory  to  obtained  the  simple  wave  solution  of  a  thin- 
walled  tube  subject  to  combined  step  loading.  Valanis  and  Wu  [39]  showed 
that  the  theory  is  able  to  predict  cyclic  creep  and  relaxation  of  metals.  Wu 
and  Lin  [40]  illustrated  how  the  strain  rate  effects  could  be  Included  in 
the  theory.  Bazant  [41],  Bazant  and  Bhat  [42],  and  Bazant  and  Krizek  [43] 
modeled  the  inelastic  properties  of  geological  materials  including  sand, 
nocks  and  concrete  using  the  theory. 

The  theory  was  attacked  by  Sandler  [44]  on  the  basis  of  a  conjecture 
that  the  theory  might  give  rise  to  numerical  instabilities  in  the  solution 
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of  wave  propagation  problems  and  also  on  its  prediction  of  unloading¬ 
reloading  behavior  which  violates  Drucker's  postulate  of  material  stability. 
To  these  arguments,  Valanls  and  Read  [45]  claimed  that  the  instability  is 
due  to  the  non-uniqueness  of  the  solution  of  posed  problem  and  not  the 
fault  of  the  endochronic  model.  They  also  claimed  that  Drucker's  postulate 
is  not  of  thermodynamic  origin  and  can  be  violated  by  the  standard 
frictional  physical  systems. 

A  more  serious  concern,  however,  was  the  openess  of  the  hysteresis 
loops.  Therefore,  the  new  endochronic  theory  was  developed  [45-47].  In  new 
theory,  a  new  intrinsic  time  was  defined  in  the  plastic  strain  3pace,  and  a 
weakly  singular  kernel  function  was  introduced.  It  was  shown  that  various 
versions  of  the  classical  plasticity  theory  are  asymptotic  cases  of  the 
endochronic  theory.  Idealized  plasticity  models  are  3hown  to  be  constitutive 
subsets  of  the  general  theory.  In  particular,  the  kinematic  model,  the 
isotropic  hardening  model,  as  well  as  their  combinations  are  derivable  from 
the  general  theory.  The  new  theory  has  been  applied  to  situations  involving 
unloading  and  cyclic  behavior  of  materials  [45,48].  Lin  and  Wu  [49]  applied 
the  theory  to  the  viscoplastic  wave  propagation  problem  of  a  thin-walled 
tube  subjected  to  impact  loading.  Valanis  and  Fan  [50]  analyzed  the  cyclic 
elastic-plastic  strain  fields  in  a  notched  plate.  Plndera  and  Herakovick 
[51]  applied  the  theory  to  model  the  response  of  unidirectional  composites 
under  off-axis  tensile  load. 

In  1978,  Valanis  proposed  an  extension  of  endochronic  theory  that  can 
apply  to  problems  of  large  plastic  deformation  [52].  He  developed  a 
constitutive  relation  for  incompressible  materials  and  used  it  to  analyze  a 
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problem  of  simple  shear.  Recently  Valanis  and  Wang  [53]  used  the 
constitutive  relation  to  analyze  the  problem  of  finite  plastic  bending. 
Further  applications  of  the  theory  to  a  set  of  special  engineering 
problems ,  including  closed  form  solutions,  will  be  presented  in  this 
dissertation. 


CHAPTER  2  -  ENDOCHRONIC  CONSTITUTIVE  RELATION  OF 


\ 


\ 


INCOMPRESSIBLE  PLASTICITY  WITH  FINITE  DEFORMATION 


An  endochronic  constitutive  relation  applicable  to  problems  of  large 
plastic  deformation  was  proposed  by  Valani3  [52]  very  much  along  the  ILnes 
of  the  endochronic  constitutive  theory  of  small  plastic  deformation  [1-2, 
46-47],  which  is  founded  on  irreversible  thermodynamics  of  internal 
variables  and  the  notion  of  intrinsic  time.  Irt  reference  [52],  the 
endochronic  constitutive  relation  of  plasticity  was  derived  for  an 
incompressible,  Isotropic,  and  isothermal  material  in  the  spatial  frame.  In 
this  chapter  the  derivation  of  the  constitutive  equation  for  large  plastic 
deformation  is  reviewed  in  the  material  frame  of  reference,  since  in  the 
next  chapter  the  theory  will  be  used  to  analyze  a  set  of  problems  which 
involve  definite  initial  configurations. 

2.1  Review  of  the  Endochronic  Constitutive  Equation 

Let  R^  be  a  material  region  with  a  surface  S^,  in  an  initial  unstressed 

state.  In  this  state  the  geometry  of  Rx  is  defined  with  respect  to  a  fixed 

Cartesian  'material'  frame,  Xa  .  The  deformed  region  Ry  with  a  surface  Sy  is 

defined  with  respect  to  Cartesian  'spatial'  frame,  Y.  .  The  deformation  of  R 

1  * 

la  defined  through  the  one-to-one  mapping  R  ♦  R  such  that 

x  y 

yl-  y^xV) 
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(2.1  .1  ) 


In  the  following  when  quantities  are  referred  to  coordinates  Xa,  their 
indices  will  be  the  lower  case  Greek  letters;  and  when  they  are  referred  to 
coordinates  ,  their  Indices  will  be  the  lower  case  english  letters. 

The  deformation  tensor  C  .  and  the  deformation  rate  tensor  d.  ,  are 

aB  ij 

defined  by 


(2.1.2) 


‘  -  <».'*>>  ha,a.  Q  .  (2.1.4) 

the  rate  of  dissipation  inequality, 

*  f*  fl  •  9  9 

Qy  -  (p,/p)  T  D  C  -  v  -  ne  2  o  ,  (2.1.5) 

as 

and  the  heat  conduction  inequality, 

-ha8,  2  0  ,  (2.1.6) 

a 


where  e  denotes  the  internal  energy  per  unit  mass, 


the  contravariant 


components  of  the  stress  tensor  in  the  material  system,  f  the  free  energy 
density  per  unit  mass,  9  the  absolute  temperature,  ha  the  heat  flux  vector, 
Q  the  rate  of  heat  supply,  Y  the  irreversible  entropy,  n  the  entropy,  p0  and 


p  the  density  of  the  medium  in  the  undeformed  and  deformed  configurations, 
respectively,  and  the  dot  denotes  the  time  derivative. 

In  the  internal  variable  formalism  an  assumption  is  made  that  the 
thermodynamic  state  of  a  body  undergoing  an  irreversible  process  can  be 
specified  by  the  current  values  of  C  and  9  as  well  as  n  internal  variables 

q  (not  necessarily  observable),  which  could  be  scalars,  vectors  or  tensors 
_r 

in  the  material  frame  or  the  spatial  frame,  depending  on  the  transformation 

laws  that  they  are  assigned  to  satisfy.  Thus  the  free  energy  ?  is  set  to  be 

a  function  of  C  ,  9,  and  q  .  Inequality  (2.1.5)  then  becomes, 

&d 


a  8  2r 


(2.1.7) 


Since  Ca0  ,  qr  ,  and  Q  are  independent  of  each  other,  the  following 

relations  must  hold  to  preserve  the  dissipation  Inequality  for  an  arbitrary 
process: 


TaS  .  2  £-  -i! 

Po  3C 


(2.1.8) 


(2.1.9) 


3?  •  .  . 

'  ?r  2  0 


(2.1.10) 


For  incompressible  materials  (  |C  j-1  ),  eq.(2.1.8)  becomes 


■  2  -g- -  P  C°» 
<30 


(2.1.11) 


where  F  is  an  arbitrary  hydrostatic  pressure. 
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The  reduced  dissipation  inequality  (2.1.10)  admits  the  internal 


constitutive  equation: 


-  —  ■  b  •  q 

3qr  -r  2 r 


(no  sum  on  r) 


(2.1.12) 


where  b  are  positive  definite  fourth  order  tensors. 

»r 

The  unique  feature  of  the  endochronic  theory  (see  Valanis  [45-47])  is 
in  the  definition  an  intrinsic  time  scale  as  the  distance  in  plastic  strain 
space  between  two  plastic  deformation  events,  and  the  stipulation  that  the 
stress  be  a  function  of  the  history  of  plastic  deformation,  measured  with 
respect  to  the  Intrinsic  time  scale.  In  the  case  of  large  plastic 
deformation  of  metals  one  can  afford  to  ignore  the  contribution  of  the 
elastic  deformation,  which  is  small,  and  thus  define  the  intrinsic  time  in 
terms  of  the  total  deformation.  Where  plastic  fluids  are  concerned,  e.g., 
metals  whose  deformation  is  so  large  that  they  have  lost  cognizance  of  their 
original  configuration  [52],  o.»e  then  defines  the  rate  of  change  of  the 
intrinsic  time  in  terms  of  the  total  deformation  rate  tensor  d^  .  More 


specifically, 


(  dt  }  “  Pijkldijdkl 


(2.1.13) 


where  should  be  an  isotropic  non-dimensional  function  of  d^  for  a 

rate  independent  material. 

The  simplest  form  of  is  a  constant  tensor 

Pljkl  ’  Pl  6lj5kl  *  p*  5ik5jl  •  (2.1.14) 

where  p,  and  p2  are  positive  non-dimensional  constants.  Eq.(2.1.13)  In 

conjunction  with  eq.(2.1.14)  yields  the  result: 
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(  &  )2 

1  dt  ‘ 


P‘diidjj  +  Pidijdij 


When  Incompressibility  applies,  dii  vanishes  and 


(2.1.15) 


(  £  )2  "  P>diJdiJ 


(2.1.16) 


An  intrinsic  time  scale  z  is  defined  to  account  for  the  hardening  or 
softening  character  of  a  material, 


2  -  d^- 
2  -  Jo  f(c,  ) 


,  (2.1.17) 

where  f(c)  is  a  positive  function  of  ?.  If  no  hardening  (or  softening)  takes 
place  f ( 0-1 . 

In  the  present  form  of  the  theory  (see  Ref. [52]),  the  free  energy 
density  is  specified  as  a  quadratic  function  of  according  to  eq.(2.1.l8): 


’  -  *  *rv,/2  c‘r)<j<j 


(2.1.18) 


where  ?0  ,  Ar,  and  C  r'  are  constants  and  q£  are  the  components  of  q  in 


the  spatial  frame. 

The  implication  of  eq.(2.1.l8)  is  that  in  the  spatial  frame  the  free 
energy  density  does  not  depend  explicitly  on  the  deformation,  a  property 
that  one  associates  with  a  "plastic  fluid",  i.e.,  a  material  that  is  so 
deformed  that  its  structure  in  the  reference  configuration  plays  no  role  in 
the  determination  of  the  Cauchy  stress.  Note  that  the  form  of  V  in 
eq.(2.1.18)  satisfies  material  objectivity  and  material  isotropy  in  the 
reference  state. 
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t 


Eq.(2.1.l8)  is  now  expressed  in  terms  of  the  covariant  components  of  q 

_r 


in  the  material  frame  and  it  is  these  that  will  play  the  role  of  independent 
variables  in  the  thermodynamic  formulation.  Thus 


*  •  ».  -Vs-  ./a 


3Y  3Y 

r  i  j  r 

where  q  na  -  — -  — ^  q.  . 

8  3X  3X8  1J 


(2.1.19) 


(2.1.19a) 


Note  that  contravariant  and  mixed  components  of  qr  In  the  material  system 


can  be  defined  by  the  relations 


a6 

3X° 

3X  8 

!(r)  ’ 

3Yi 

55J ' 

a 

3X° 

(r)B 

“  3Yi 

3X8 

(r )  8 

3X8 

1  a 

3Xa 

3Y 

(2.1.19b) 


(2.1.19c) 


(2.1 .I9d) 


The  internal  constitutive  equation  (2.1.12)  is  now  written  in  the 
specific  form 


h(r)  !(r)  3?  . 

b  qc6  +  7^6  "  ° 

3q(r) 


(2.1.20) 


(r) 


where  b  are  scalars  and  the  roof  denotes  the  derivative  with  respect  to  z. 


It  implies  that  the  endochronic  rate  of  change  of  the  covariant  components 

( r ) 

of  9  should  be  proportional  to  the  covariant  components  of  the  internal 


1 


3?  (r*) 

stress  tensor  — -  ,  in  which  case  b  is  of  the  form: 
a0  ■ 

q(r) 

h(r)  h(r) 
baBY<5  b  5aY  5S6 


(2.1.20a) 


The  solution  of  eq.(2.1.20)  in  conjunction  with  the  form  of  f  in 

eq. (2.1 .19)  and  the  initial  condition  on  q  .  (q  L«  <5  „),  dictated  by  the 

aB  aB'  aB 

fact  that  the  Cauchy  stress  in  the  reference  state  can  be  at  most  a 
hydrostatic  pressure,  is  given  by 

(r)  A(r)  „  A(r)  *  -a_(z-z’):  ... 


F7caa*77T^e~V  c<.8(z'>“'  •  (2-’-21) 


where  a 


r  b(r) 

On  the  other  hand  eq. (2.1.11)  in  conjunction  with  eq.(2.1.l9)  gives 


aB  -aB  ,,,  (r)  aB  ,  r (r)  (r)a  (r) 8’ B. 
T  --PC  2(4  V)-  C  q  s,q  ) 


(2.1.22) 


Ta8  P  Ca8  2(A  qaB  C  qaB'qB  5  *  (2’K23) 

Substitution  of  eq.(2.1.21)  into  eq.(2.1.23)  gives  the  following 
relation  for  the  covariant  components  of  the  stress 


T»«-  -pe.e*  - 


’cU  ’l  ! 0  izx&zt 


.  (2.1.24) 


It  was  found  from  an  analysis  of  the  problems  of  simple  extension  and 
simple  shearing  that  the  contribution  to  the  stress  of  the  second  integral 

relative  to  that  of  the  first  integral  is  of  the  order  of  1/a2  where  a  is  of 
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the  order  of  one  hundred,  this  being  the  ratio  of  the  elastic  modulus  to  the 
maximum  stress  (in  tension  or  in  shear).  In  the  subsequent  analysis  the 
double  integral  was  ignored  with  the  result  that  the  covariant  components  of 
stress  in  the  material  frame  of  reference  are  given  in  terms  of  a  linear 
history  integral  of  the  right  Cauchy-Green  tensor  Cag.  Thus 


T  «  -P  C  0+  J^G(z-z’)  C  Lz')  dz' 
a&  aB  0  aB 


(2.1.25) 


where 


.<r)  -  , 

G(s)  “  2  l  -yyy  e  ar 


(2.1.25a) 


For  an  extensive  discussion  see  Ref. [52]. 


2.2  Description  of  Constitutive  Equation  in  Curvilinear  Systems 

s 

In  order  to  describe  the  deformation  in  a  curvilinear  system,  we  set  4> 
and  to  be  curvilinear  coordinates  for  the  material  system  Xa  and  the 
spatial  system  ,  respectively,  as  shown  in  Fig.l. 

For  tne  material  system,  we  have  the  relation, 


♦  8-  *8(Xa) 


(2.2.1) 


The  corresponding  covariant  and  contravariant  metric  tensors  and  G  are 
given  by  eq.  (2.2.2) : 


5Xk  3Xk 
aS"  3*a  3$8 


(2.2.2a) 


Ga  •  ( T / { G | )  (cofactor  of  Gog) 


(2.2.2b) 
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where  j  G {  i3  the  determinant  of  Gag. 

The  corresponding  relations  for  the  spatial  frame  are 


v  vv 


3Tl 

®kl“  "30"  30” 


(2.2.3) 


(2.2.4a) 


and 

gkl«  ( 1  /  | g|  )  (cofactor  of  g^)  (2.2.4b) 

The  deformation  map  (2.1.1)  can  be  then  written  through  eq.(2.2.1)  and 
(2.2.3)  in  terms  of  one  of  the  following  alternatives,  depending  on  the 
problem  at  hand: 

Yi-  f1(eB,t)  ,  (2.2.5) 

9j*  9j(xa,t)  (2.2.6) 

and 

V  ej(4S,t)  (2.2.7) 

The  Cauchy-Green  deformation  tensor  C  and  its  curvilinear  counterpart 


C  are  given  below: 
uv 


c  !!i 

olB  gfclgXa  3XS 


(2.2.8) 


or  in  matrix  form 

cci  -  cnT  [g]  [f] 


(2.2.3a) 


where  the  deformation  gradient  matrix  [F]  is  given  by  eq. (2.2.9). 


F  -  — - 

ka  3Xa 


(2.2.9) 


also. 


r  .  .  !fi 

“v  kl  3*“  3«v 


(2.2.10) 


or 


[C]  -  Cf V  [g]  CF] 


(2.2.10a) 


where 


}  .  Sl 
ku 


(2.2.11) 


We  now  prove  the  following  relatione,  as  these  will  be  of  Interest  in 
what  follows: 

(1)  The  incompressibility  condition  is  given  by 


|C  I  -  |G.  .  j 
1  uv1  1  kl* 

Proof:  From  the  tensor  transformation. 


(2.2.12) 


c  .c  £! »! 

“v  a8  3," 


(2.2.13) 


Then 


•  I  <S>TI-|CJ-|^I-  l»TH»l 


(2.2.14) 


where  jC^  [  -  1  is  unity 


and 


M 


3XJ 

3$x 


(2.2.15) 


The  proof  is  completld  by  noting  eq.(2.2.2a). 


(2)  The  deformation-rate  tensor  d^  is  also  given  by 


d  i/p  m!  c 
dij  /2  SYj  CaS 

proof:  The  rate  of  change  of  C  .is 

as 


(2.2.16) 


3v  3Y,  3Y  3v. 
£  _  1_  i_  +  i_  i_ 

aB  3Xa  3X8  3Xa  3XS 


3v  3v,  3Y  3Y. 

-A  +  -J-  )  _1  -J- 

3*j  3Yi  3Xa  3X8 


3Y  3Y 

m  2  d  — *•  — 

3X°  3XS 


(2.2.17) 


from  which  eq.(2.2.l6)  follows  easily. 


(3)  The  intrinsic  time  scale  can  be  obtained  by  means  of  the  following 
relation,  derived  by  direct  substitution  of  (2.2.16)  into  eq.(2.l.l6). 

2£  i* 


(  )*  -  (p,  c;1  c  .  c;i  cnY)/4 

dt  *  Ya  a8  8n  nY 


(2.2.18) 


or 


(  dc  )l  -  (patr([D]))/4 


(2.2.19) 


where 


[D]  -  [C]"1  dCG]  CC]'1  dCC] 


(2.2.20) 


In  the  same  manner,  we  also  can  obtain  d;  from  eq.(2.2.12) 


(  dc  )*  -  (patr([D]))/4 


(2.2.21 ) 


where 


CD]  -  Ccf 1  dCG]  CC]"1  dCC] 


(2.2.22) 


18  - 


Coordinate  systeaa 


20 


CHAPTER  3  “  CLOSED  FORM  SOLUTIONS 


To  illustrate  the  application  of  the  theory  to  the  domain  of  finite 
plastic  deformation,  in  this  chapter  the  theory  will  be  used  to  analyze  a 
set  of  special  and  yet  practical  problems.  These  problems  are: 

(1)  large  plastic  uniform  extension  of  a  cuboid; 

(2)  large  plastic  shear  flow  in  a  circular  tube; 

(3)  large  plastic  torsion  of  a  circular  bar; 

(4)  large  plastic  bending  of  a  block; 

(5)  large  plastic  expansion  of  a  sphere. 

Closed  form  solutions  for  all  problems  are  obtained.  To  our  knowledge, 
this  is  the  first  time  in  the  field  that  this  set  of  problem  are  solved  in 
closed  form  solutions.  In  the  course  of  obtaining  the  solutions  the  semi¬ 
inverse  method  is  used.  The  deformation  field  is  derived  from  the  condition 
of  incompressibility  to  within  a  set  of  unknown  parameters,  which  are  then 
determined  by  satisfying  the  equilibrium  and  boundary  conditions. 

In  order  to  find  the  complete  solutions  for  these  problems,  we  need  to 
know  the  kernel  function  G(z)  in  the  constitutive  equations  (2.1.25), 
(2.2.24),  and  (2.2.25).  In  the  development  of  the  endochronlc  theory  of 
plasticity  in  the  region  of  infinitesimal  deformation,  the  kernel  function 
was  discussed  thoroughly.  There,  G(z)  was  required  to  be  weakly  singular  at 

the  origin  but  integrable  in  the  domain  0SzS«,  i.e.,  G(0)-«  and  / oG(z)dz<®. 

This  requirement  leads  to  the  following  consequences:  (1)  it  gives  rise  to 
closed  hysteresis  loops,  (2)  it  ensures  initially  elastic  unloading,  i.e., 
Initially  zero  rate  of  dissipation  upon  unloading,  but  (3)  it  implies  an 


-  21 


ok 

• 


infinitesimally  small  elastic  domain.  Two  explicit  algebraic  forms  were 
discussed  at  great  length  in,  [45].  They  are 


(i) 


P  0 

G(  2 )  -  — 


I  0p.'V 


(3.D 


where  0<<n<1 ,  p„>0,  Gr>0,  and  for  all  r,  and 

09 

(li)  G( 2 )  -  I  G  e~arZ  (3.2) 

i  r 

with  the  conditions  that  G  and  a  be  positive  for  all  r, 

r  r 

A 

I  G  -  •  ,  (3.3a) 

i  r 

and 

a 

l  G  /a  <  -  (3.3a) 

1  1  r 

The  determination  of  the  material  functions  is  very  important  in  order 
that  the  plasticity  theory  can  be  developed  and  applied.  It  is  a  task 
involving  a  considerable  amount  of  experimental  investigation  and 
theoretical  analysis.  This  is  not  a  object  of  the  present  dissertation.  In  a 
recent  work  [54],  the  kernel  function  was  determined  by  a  strain  controlled 
cyclic  test.  The  kernel  function  G(t)  was  given  by  the  slope  of  the  cyclic 
plastic  strain-stress  curve  for  the  shear  test  at  a  steady  state.  The 
hardening  function  f(c)  of  eq .(2.1.17)  was  determined  from  reversals  of  the 
cyclic  shear  test.  Fan  [55]  determined  by  an  approximate  method  the  material 
functions  which  were  convenient  for  engineering  application. 
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In  this  work,  to  demonstrate  the  application  of  the  theory,  eq.(3.i)  or 
eq.(3.2)  are  not  applied  to  a  specific  material  in  detail.  Instead,  just  for 
a  convenience  of  calculation,  an  approximate  form 

G(z)-G,e‘aZ  (3.4) 

will  be  taken  in  problems  of  uniform  extension,  bending,  and  expansion, 
while 

G(  z)-pl)/zUi  (3.5) 

in  the  problems  of  shear  flow  and  torsion.  These  kernel  functions  represent 
simple  and  yet  realistic  cases  [52,53].  For  numerical  calculation,  a  will  be 
set  to  200  (this  number  is  characteristic  of  pure  aluminum  with  a  tensile 

modulus  of  107  lb/in2  and  an  ultimate  stress  of  50*1 03lb/in2 ,  a  being  their 
ratio.)  and  w  to  0.86  (from  [48]  for  normalized  mild  steel). 

3.1.  Large  Uniform  Plastic  Extension 

Uniform  extension  of  a  unit  cube  is  considered.  Taking  Cartesian 

coordinates  x°  and  y^  (which  parallel  the  sides  of  the  undeformed  and 

deformed  bodies)  as  the  material  and  the  spatial  coordinates,  respectively, 
the  deformation  field  can  be  expressed  by  the  following  relations: 

y,  -  l,(t)x‘,  y2  -  X2 ( t) x2 ,  y,  -  x,(t)xJ  (3-1.1) 

where  X,,  X2,  and  X,  are  the  stretches  in  the  xl,  x2,  and  x1  directions 
respectively,  X^>1  corresponding  to  extension  and  0<X^<1  to  compression. 
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For  a  particular  case  of  simple  extension  under  a  force  parallel  to  the  xl 
direction  of  the  cube,  as  shown  in  Fig. 2,  we  have  X,>1  and  0<X2-Xj<1 . 


The  Cauchy-Green  deformation  tensor  is 
X?  0  0 


[C] 


OOXl 


(3.1.2) 


1/2 


For  Incompressible  material  | C | —  X*  XJ  -  1,  i.e.  X,-  (1/X,)  .  Hence  the 

deformation  for  an  incompressible  cube  becomes 

y,  -  X ,xl  ,  y2  -  (1/X,)1/2xa  ,  y,  -  (1/X,)1/2x*  .  (3-1. 3) 

The  Cauchy-Green  deformation  tensor  now  becomes 


[C] 


X?  0 


0  1/X,  0 


0  1/X, 


(3.1 .«) 


Then  the  inverse  and  the  increment  of  Cauchy-Green  tensor  are  given  by. 
1/X?  0  0 


[C] 


-1 


X  i 


0  X, 


(3.1.5) 


and 


With  direct  substitution  of  eq.(3-1.5)  and  eq.(3.1.6)  into  eq.(2.2.20),  the 
intrinsic  time  is  obtained  as 

(dc)1-  3pa/2  (dX,/X,)a  .  (3-1.7) 

Setting  pa  equal  to  2/3  (the  actual  value  of  p2  is  immaterial  for  all 

problems  in  this  paper),  (dO*»  (dX,/X,)*.  For  monotonlc  extension, 

dc-dX,/X,  .  Upon  Integration,  c-ln(X , )+c,  where  c  is  an  arbitrary  constant. 

The  condition  c-0  for  Xt-1  requires  that  c-0  and  c-lnU,).  Neglecting 

hardening  (or  softening)  effects  in  the  material,  f(c)-1  and  dz  and  z  are 
obtained  as  shown  in  eq. (3.1.8). 

dz-dXj/X,  ,  z-ln(X,)  .  (3.1.8) 

With  the  help  of  eq -(3.1 -3)  [C]  is  given  below, 

2X  *  0  0 

[C]  -  0  -X?1  0  (3.1.9) 

o  o  -xT1 

‘  . 

-CXZ 

Setting  G-Gae  ,  the  stresses  are  obtained  upon  use  of  the  constitutive 
equation  (2.1.25)  as  follows. 


T,,-  X?  [-P 


(3-1 -10a) 


^VJ2  nt+O 

*  -ozzt  (1 '  ux‘  11 
G« 

Tl2-  Tn  -  i/x,  c-p  -  JZ?T)  0  ■  1/x>  )3 

and 

‘rU-T2J-Tlj"0 


(3.1.10b) 


(3.1.10c) 


The  components  of  the  stresses  in  spatial  coordinates  are  given  by 
T»i-  -P*  j^2)  (1  -  1/X?"2)  ,  (3. i.lla) 

T22-  T,,  -  -P  — (1  -  I/Xf1)  (3.1.11b) 

and 

Ttl-  T,,  -  T»,  -  0  .  (3.1.110 

In  the  absence  of  body  forces,  equilibrium  in  the  deformed  body  with  the 
stresses  in  eq .(3.1.11)  gives  3p/3y,-  3P/3y2-  3p/3y,»  0  ,  i.e.,  P  is  a 

constant.  This  constant  is  determined  by  satisfying  the  boundary  conditions. 
In  the  present  case,  T22  and  T,,  are  zero  at  the  edges  in  y2  and  y3 

directions.  Hence, 

G« 

p  -  -  (i  -  i/*r  )  .  (3.i.i2) 

Substituting  eq.(3.1.12)  into  eq. (3.1.11),  the  only  no-vanishing  stress  T^ 
(normalized  by  G0)  is  given  below. 


tii/Go  .  0  -  ^/xTz)  *  -±r  (1  -  . 

while  the  resultant  force  in  direction  is  given  by 


(3.1.13) 


3a 


F/O.  -T..J5/0.  '  .  (3.1.1*) 
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When  X,  becomes  very  large  ,  T,,  approaches  a  constant  given  by 


3aGe 


Tt: 


(3.1.12) 


1  (a+2)  (a-1 ) 

However,  since  the  cross  sectional  area  of  the  cuboid  decreases  with  Xi ,  the 


resultant  force  will  decrease  after  a  critical  value  of  X*  and  will  go  to 
zero  as  X,  goes  to  infinity.  The  variation  of  Tn  and  F  with  Xj  are  3hown  in 


Fig. 3  and 


3.2.  Large  plastic  flow  in  a  circular  tube  (pipe  flow) 

Consider  a  case  of  pipe  flow  in  which  each  point  in  the  material  moves 
parallel  to  the  axis  of  the  pipe  through  a  distance  f(r)  depending  only  upon 
the  radial  position  of  the  point. 


Taking  the  cylindrical  polar  coordinates  $  (R,0,Z)  in  the  undeformed 


body  and  ©^(r.e.w)  in  the  deformed  body,  the  deformation  field  is  expressed 


as  follows: 


r»R,  9-0,  w«-K(t)f(r)+Z  .  (3.2.1) 

where  K  is  a  positive  function  and  the  deformation  corresponding  to  positive 
stresses  is  shown  in  Fig. 5. 

The  metric  tensor  in  the  material  system  is 


1  0  0 


3 


0  R*  0 
0  0  1 


lGJ-R2 


(3.2.2) 
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and  in  spatial  system 


[g  ]  - 

U8klJ 


1  0  0 
0  r*  0 
0  0  1 


|Skil“r2  *  0.2.3) 


The  deformation  gradient  is 

1  0  0 


CF]  - 


1 


-Kf’  0  1 


(3.2.4) 


where  f'-df/dr.  The  Cauchy-Green  deformation  tensor,  derived  by  use  of 
eq. (2.2.10a) ,  is 


\  1 +K Jf ' 1  0  -Kf'  I 


0 


-Kf’ 


r1  0 
0  1 


|C  | *rl  .  (3.2.5) 

1  uv 1 


The  incompressibility  condition  is  satisfied  automatically  since 

|C  I-Ig,  ,1.  The  inverse  and  the  Increment  of  [C  ]  are 
'  uv1  1  kl '  uv 


1  0  Kf’ 


0  1/r*  0 


Kf’  0  1 ♦K2f ' * 


and 


2Kf'*dK  0  -f’dK  1 


dCC  ]  - 


0  0  0 


(3.2.6) 


(3.2.7) 


The  Intrinsic  time  can  be  found  from  eq. (2.2.22)  to  be 


(dO2-  (f'dK)* 


(3.2.8) 


*g£r,  where  p*  has  been  set  to  2.  For  monotonic  shearing,  d?-f'd K.  After  an 


S  integration  and  imposition  of  the  initial  condition,  c-0,  at  K-0,  it  follows 

■4 

that  c-f'K.  Upon  neglecting  the  effect  of  hardening  (or  softening)  of  the 


'  Hi  material,  d z  and  z  are 

# 


dz-f’dK  and  z-f’K 


(3.2.9) 


1 


^  Now  The  Cauchy-Green  deformation  tensor  is  rewritten  with  the  help  of 
rjf  eq. (3-2.9)  in  the  form 


1+z2  0  -z 


[C  ]  -  0  r2  0 

uv 


-z  0  1 


and  d[C]/dz  is  obtained  in  the  form 


2z  0  -1 


(3.2.10)) 


dCC]/dz 


0  0  0 


(3-2.11 ) 


[  -1  0  0  J 

Setting  G(z)  -  p0/z°,  the  stresses  are  obtained  by  use  of  the  constitutive 


equation  (2.2.25)  in  the  form 


tRR— P(1>z2)  +  2p0/(1-a)/(2-a)  z 


(2-a) 


(3.2.12a) 


tP7-Pz-‘p,/(  1-a)  z 


( l-o) 


(3.2.12b) 


T99"”Pr’J  ’  TZZ^P  ’  tR0“  T0Z“  0 


(3.2.12c) 
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The  components  of  stress  in  the  spatial  curvilinear  system  T  are  given  by 


frr-  -P-2p./(2-.)  .  fM-Pr*  .  f^-P  .  (3_2_,3) 

fru-  -p»/(1->a)  i(t  “)  ana  Tpe-  Teu-  0 

while  the  physical  components  of  stress  in  the  spatial  system  are  given  by 


°rr"  "P+g‘  ’  °89"  °ww“  "P  ' 


a  ■  g, 
rw 


and 


o  •  a.  •  0 
re  8w 


where 


g,-  -2p#/(2-o)  z(2  °I  -2p0/(2-o)  (f'K)(2‘a) 


(3.2.14) 


(3.2.14a) 


and 

ga-  -p./O-o)  z(1^a)-  -p,/(1-a)(f’K)(1*o)  (3.2.14b) 

In  the  absence  of  body  forces  the  equations  of  equilibrium  along  with  the 
stresses  of  eqs.(3.2.l4)  give  the  relations, 

-3P/3r  +  3g,/3r  +  gt/r  -  0  ,  (3.2.15a) 

-3P/3w  ♦  3ga/3r  +  g2/r  -  0  (3.2.15b) 

and 

-3P/39  -  0  (3.2.15c) 

The  last  equation  requires  that  P  be  a  function  of  r  and  w  only.  Since  g, 


depends  on  r  only,  equation  (3.2.15a)  requires  32P/3r3w-0.  Thus, 
P  -  C„w  ♦  Y(r)  ,  (3.2.16) 

where  C0  is  a  constant  of  integration.  From  equation  (3.2.15b), 

g*  -C ,/r  ♦  C,r/2  ,  (3.2.17) 
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where  Cx  Is  another  constant  of  integration.  To  ensure  a  finite  solution  at 


1 


Upon  satisfying  the  boundary  conditions  (i)  and  (ii)  in  eq.  (3.2.23)  the 
following  equations  are  obtained: 


-f  »  -Dtra 


^»p,(3-2tt)(1"tt)  j(2-qt)/(1-g)a(4- 3tt)/(1-a) 


(2-a)2  (4-3a) 


(3.2.24) 


.  „  ,  „  .  ,  4p#(3'2o)(1‘a),“C9(1'a)i(2-a)/(1-a).(4-3o)/(1-a) 

)--DTra2-C0L.ira2+—  —  —  •{ — = - }  a 


(2-a)2  (4-3a) 


Solving  the  simultaneous  equations  (3.2.24),  The  constants  C„  and  D  are 


obtained  as 


C„  -  -F/  (Lira2 ) 


(3.2.25) 


4p o (3— 2a )  (1  —a)  — C0(1— a)  /.  w. »  ( o— 

-  F/(ira*)  ♦  -  {  -  j(2-a)(!-a)r(2-a)(T-a)  # 


(2-a) 2 (4- 3a) 


Finally,  the  following  non-vanishing  stresses  are  obtained  and  are  given  by 
eq. (3.2.26) 


°rw/p°"  •0(r/L)/2 


(3.2.26a) 


°ww/p#“  a08/po"  8(w/L-1  )-2(3"'2a)/(2-a)*[(  1-a)/2*8] 


(2-a)/(1-a) 


•{2(1-a)/(4-3a)(a/L)(2"a)/(1"a)-(r/L)(2"o)/(1"a)} 


°rr/p,“  9ww/P«-'2/(2-a)[(1-a)/2-3](2"a)/(1'a)(r/L)(2"a)/0'a) 


(3.2.26b) 


where 


(3.2.26c) 


8  -  F/(irp„a2) 


(3.2.26) 
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The  displacement  in  the  w  direction,  u,-  w-Z  -  -Kf(r),  can  be  obtained 

by  an  integration  of  the  expression  z-Kf’(r)  in  eq.(3.2.20)  and  is  given  by 

Kf-CF(t-a)/(2LiraJp0]1/(1‘a)[(1-o)/(2-a)]r(2'°)/(1"o)^  C,  .  (3-2.27) 

Imposing  the  boundary  condition  (Hi)  in  eq. (3.2.23)  gives  the  integration 
constant  C,, 

C,-  -[F(1-o)/(2Lira*p0]1/(1_a)[(1-a)/(2-a)]a(2'a)/(1"a)  (3-2.28) 


Hence  the  displacement  u, ,  normalized  by  a,  the  radius  of  pipe,  becomes 


C(1-a)/2*8]1/(7"o){(r/a)(r/L)1/(1"a-(a/L)1/(1"a)}  . 


(3-2.29) 


It  is  noticed  that  us  reaches  a  maximum  value  at  r-0,  which  is 


u,(0)/a  -  4!”r  Ca(1-a)/2L]1/(1‘a)S1/(1‘a) 


(3-2.30) 


The  above  solutions  can  be  used  to  analyze  the  process  of  metal 
extrusion  in  a  cylindrical  die.  The  condition  to  move  a  bar  of  length  L  and 
radius  a  is  given  by  <Jrw“  (<Jpw)c  at  r"a-  (<Jrw)c  tlie  factor  representing 

the  surface  condition  between  the  metal  and  the  die.  When  F  increases  until 

F  ,  o  reaches  the  critical  value  (<j  )  at  r-a.  The  metal  bar  starts  to 
v-  rw  rw  c 

oove  in  the  die  with  the  shape  of  u,  corresponding  to  F  . 


(3.2.31a) 


C(1-o)/2.80]1/(1'a){(r/a)(r/L)1/(1'a^(a/L)1/(1'a)} 


where 


(3.2.31b) 


8C-  Fc/(TTp,a*) 
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With  a-  0.86  a  numerical  example  la  given  aa  followa.  Fig. 6  shows  the 
deformation  profile  for  varioua  parameter  g.  Fig. 7  givea  the  relation 
between  deformation  u,  and  the  parameter  8  ,  i.e.,  the  applied  force,  at 

r-0.  o  ia  a  linear  diatributlon  over  the  croaa  aection  of  the  pipe,  a  -0 
*  **  rw 

at  r-0  and  a  -  8a/ 2L  at  r-a.  The  diatribution  of  o  and  <j  over  the  croaa 

4  w  ww  w 

aectlona  w-0,  and  w/L-0.5,  are  given  in  Fig. 8  for  the  caae  of  a/L-0.2,  and 
6-70. 

3.3.  Large  plaatic  toraion  of  circular  bar 

Conaider  a  uniform  aolid  circular  bar  of  radius  a.  The  one  end  of  the 
bar  ia  fixed.  The  other  end  la  subjected  to  an  angle  of  twiat  due  to  applied 
torque  T.  The  bar  ia  also  assumed  to  be  constrained  axially,  thus  allowing 
the  possible  development  of  an  axial  force  F.  The  deformation  field  ia  the 


following: 


r-R  , 

e-e+k(t)z  (3.3.1) 

and  w-Z  , 

where  iR,9,Z)  and  (r,8,w)  are  the  cylindrical  coordinates  for  the  material 
and  spatial  system,  respectively.  Eq .(3.3.1)  implies  that  the  planes 
perpendicular  to  the  axis  of  the  bar  are  rotated  in  their  own  planes  through 
an  angle  proportional  to  the  distance  of  the  plane  from  the  fixed  end,  k  ia 
the  twiat  per  unit  length. 

The  metric  tensor  in  the  material  system  ia 

1  00' 


(3.3.2) 


-  3«  - 


'  1  ■  ' II  |(i|jip 

\ 
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i:  The  intrinsic  time  can  be  obtained  as 

V? 

(dc)2-  (rk)2  ,  (3-3.8) 

where  P2  has  been  set  to  2.  Under  the  condition  of  monotonio  torsion, 


dc-t*dk.  After  integration  and  consideration  of  a  zero  strain  initial 
condition,  c-rk.  Neglecting  hardening  and  softening  effects,  dz  and  z  are 
given  below, 


dz-rdk, 


z-rk 


(3.3.9) 


Hence,  the  Cauchy-Green  deformation  tensor  can  be  written  as 

1  0  0 


CC  ] 

yv 


Then, 


[C  ]  - 
yv 


0  r2  rz 

0  rz  z2+1 

0  0  0 

0  0  r 

Or  2z 


(3-3.10) 


(3.3.1D 


Again  taking  G(z)  -  p0/z  ,  the  following  stresses  are  obtained  after 

substitution  of  eq .  ( 3. 3'.  1 0)  and  eq .  ( 3 . 3 . 1 1 )  into  eq.(2.2.25)  and  then 
integrating  the  resultant  equations, 
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trr"  “p  ’  Tee“  “Pr*’'  tR0*  trz“  0  * 


tZ2-p(z2-i)- 


2p, 

( 1-a) (2-a) 


t 

£  and 


Tgz  -  -Prz+[p0/(1-a)]  t*2 


(2-a) 


( 1-a) 


(3.3.12) 


The  stresses  in  the  spatial  curvilinear  system  T  can  be  obtained  by  tensor 
transformation  as  follows: 

(2-a 


Trr"  _P  ’  T00*"Pra  ’  Tww-P-C2p0/(2-a)]zv 

f9w-  Cpa/(1-o)]-r.2(1'o)  and  TpQ-  T^-  0 


(3.3.13) 


The  physical  components  of  the  stresses  in  the  spatial  system  are 
°rr“  °eQm  ~?  ’  ffww“  -1 -[2p0/(2-a)]  z'2  o) , 

l 

°0w  “  EPo^1***)]  z(1  a)and  ar0-  apw-  0  ( 3 - 3 - 1  ^ ) 

When  body  forces  are  zero,  the  equilibrium  equations  with  the  stresses  in 
eq.G^.I1*)  give  3p/3r-0,  3p/30-O,  and  3p/3z-0.  Thi3  states  that  P  is  a 
constant.  If  the  surface  of  the  cylinder  r-a  is  to  be  free  of  tractions, 
°rr"°  there.  Hence  P  must  be  equal  to  zero.  Consequently,  the  following  non- 

vanishing  stresses  yield 

(1-a) 


°0w  *  Cp0/(1-a)](rk] 


(3.3.15a) 


®ww“  C-2p,/(2-o)](rk) 


(2-a) 


(3.3.15b) 

Clearly  there  is  no  any  integration  constant  available,  the  forces  at  the 
ends  can  no  longer  be  controlled  for  the  prescribed  deformation  field.  The 
resultant  force  F  which  is  needed  to  maintain  the  length  of  the  bar  during 
the  deformation  can  be  calculated  as 
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(3.3.16) 


F/(p«a*)  -  (;^|aww|2irrdr)/(p0aJ)  -  (2-cO(4-a)  (ka)(2  a) 
while  the  resultant  torque  as 
T/(Poa* )  -  (/^2ir00wradr)/(p 9a1 )  -  (ka)(1~a)  .  (3-3.17) 

Solving  for  k  from  (3.3.17),  the  twist  per  unit  length  is  obtained  and  given 
by  eq. (3-3.18), 

r~,,  .ii  ( 1_a) (4-a)  1 1/(  1-a)  , 

k  -  1/a  [T/(poaJ)  - ^ - J  (3.3.13) 

The  total  twist  angle  at  any  station  w  can  be  obtained  as 

A9  -  0-e  -  kz  -  (w/a)[T/(p0a*)  l-1.~?^(4~a),]1/(  1~ct)  .  (3.3.19) 

The  numerical  results  with  a-0.86  are  given  in  Fig. 10  -  Fig.  13-  Fig.  10  and 
Fig. 11  show  the  variations  of  the  torque  and  axial  force  vs.  the  twist, 
respectively.  The  distributions  of  0  .  and  a  over  the  cross  section  are 

W0  ww 

given  in  Fig.  12  and  Fig. 13.  It  is  interesting  to  observe  that  for  relative 
small  k  ( k < < 1 )  the  torque  increases  more  rapidly  than  the  axial  force, 
however  for  relative  large  k  (k>1)  the  axial  force  Increases  much  faster 
than  the  torque.  For  the  elastic  case  when  o-O,  in  eq. (3-3.17)  and 
eq. (3.3.16) 

T  -  (irp#aV2)k  ,  F  -  (irp,aV2)k*  (3-3.20) 

For  infinitesimal  deformation,  k<<1  ,  F  vanishes.  The  results  agree  with 
those  from  the  theory  of  elementary  mechanics  of  materials. 


3.4.  Large  Plastic  Bending  of  a  Block 

Let  a  block,  bounded  by  the  planes  X}-  a,,  X|-  a2 ;  and  X1-  ±b;  X*-  ±c 

In  the  undeformed  state,  be  deformed  symmetrically  with  respect  to  X1  axis 
into  a  portion  of  a  cylinder  as  shown  in  Fig. 14.  Letting  0j(r,0,w)  be 

coordinates  in  deformed  state,  the  deformation  map  can  be  stipulated  by 

r  -  r(X‘ ,t)  , 


0  -  0(Xa,t) 


(3.4.1) 


and  w  -  w(X*,t) 

It  will  be  shown  presently  that  the  condition  of  incompressibility  restricts 

the  deformation  given  by  eq.(3.4.1)  to  a  form  where  a  plane  of  constant  X1 

deforms  to  a  surface  of  constant  r;  a  plane  of  constant  X2  to  a  plane  of 

constant  8;  and  a  plane  of  constant  X*  to  a  plane  of  constant  w. 

The  metric  tensors  in  the  spatial  system  are 


1  0  0 


1  0  0 


^l"-1 


0  r2  0 


[g  3-0  1/r2  0 


(3-4.2) 


0  0  1 


0  0  1 


The  deformation  gradient  matrix  is 


dr/dX' 


d0/dX* 


dw/dX1 


-  39  - 


0  -  c,Xl  ♦  A,(t)  (3.4.6) 

^  w  -  c,X*  ♦  A,(t) 

where  A,,  A,,  and  A,  are  the  constants  of  the  integrals. 

Considering  the  plane  strain  condition  (o,-1,  a,«0)  and  excluding  rigid 
body  rotation  (A,-0),  the  deformation  field  becomes 


r*/2  -  R  X1*  A(t)  , 

0  -  (l/R)  X*  (3.«.7) 

and  w  -  X1 
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y**- 

is* 


where  the  facts  that  c*  is  the  curvature  and  c1  the  radius  of  curvature  of 

the  deformed  cylinder  have  been  used.  The  Cauchy-Green  tensor  and  its 
inverse  now  become 


(3.4.8) 


1/Q* 

0 

0  ' 

[c  f1  - 

'  Q* 

0 

0  ' 

[C]  - 

0 

Q* 

0 

and 

0 

1/Qa 

0 

0 

0 

1 

0 

t 

0 

1 

where  Q  -  r/R.  The  increment  of  [C]  is 

0  0 

d[C]  -  0  2QdQ  0 


(3.4.9) 


0  0  0 
The  intrinsic  time  can  be  obtained  from  eq.(2.2.15)  and  eq.(2.2.l6)  whereby 


(dc)*  -  C(1/Q)dQ]*  , 

where  p2  in  eq.(2.2.15)  has  been  set  equal  to  1/2. 


(3.4.10) 


For  monotonic  bending  dc  -  (l/Q)dQ  and  then  c  *  ln(Q)*c  after  integration. 
Since  c-0,  for  Q-1,  hence  c-0.  Neglecting  hardening  or  softening  effects, 
f-1  and  dz  and  z  can  be  derived  as  shown  in  eq.(3.4.11), 

dz  -  ( 1/Q)dQ  ,  z  -  ln( Q)  .  (3.4.11) 

With  the  help  of  (3.4.11)  C  is  given  below: 


CC]  -  Q 


-2Q 

0 

0 


-3 


0 

2Q 

0 


0 

0 

0 


(3.4.12) 
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taking  G(z)  -  Cte  02 ,  the  stresses  are  obtained,  upon  use  of  the 
constitutive  eq.(2.1.25), 


T,,-  (1/Q*)C“P-(2G0/(a-2))(1-l/Q(a-l))] 

TJ2-  Qa[-P+(2G0/(a*2))(1-1/Q(o+2))]  ,  (3-4.13) 

Tj j*  —  P ,  and  Tj j*  T  2  2  •  T j j*0  - 

The  components  of  stress  in  the  spatial  curvilinear  system  are  given  by 
eq. (3*4.14), 


T  -  -P  -(2G0/(a-2))(1-1/Q(tX_2))  , 

rr 


t9Q-  ra[-P+(2G0/(a-*-2))(1->1/Q(a+2))]  , 


T  -  -P  ,  and  T  -  T  -  T  -  0 
ww  re  rw  we 


(3-4.14) 


while  the  physical  components  of  stress  in  the  spatial  system  are  given  by 

eq. (3-4.15), 


°rr“  -P-(2(V(a-2)(1-1/Q(o'2)) 


a09-  -P*(2G,/(a+2)(1-1/Q(a+2)) 


a  -  -P  and  a  -  a  -  0.  -  0 
ww  re  wr  9w 


(3-4.15) 


where  P  will  be  determined  from  the  equilibrium  equations  and  the  boundary 
conditions. 

When  body  forces  are  absent,  two  of  the  three  equations  of  equilibrium, 
in  the  e  and  w  directions,  are  satisfied  identically  if  the  stresses  in 
(3-4.15)  are  functions  of  r  only.  The  remaining  equilibrium  condition  is 
given  by  eq. (3-4.16): 
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(3.4.16) 


rr 


+ 


30 


0 


Substituting  the  expressions  for  the  stresses  in  eq.(3.4.15)  into 
eq. (3*4.1 6)  and  integrating,  the  following  relation  for  P  is  obtained, 

20 

P  . - {n(R/r)(““2)}  +2G0{-2a/(a4-4)ln(r/c)-1/(a-2)2(R/r)(a"2)  } 

cl- 2 

-1/(a+2)*(R/r)(a+2)}  ,  (3-4.17) 

where  c  is  the  constant  of  Integration. 

Normalized  by  2G0,  the  non-vanishing  stresses  in  (3.4.15)  with  the  help  of 


eq. (3-4.17)  become 

a  /2G0-  2a/(aI-4)ln(r/c)  +  1/(a-*2)a(R/r)(°"‘2)H/(a+2)2(R/r)(a+2).  (3.4. 18a) 
rr 

<jqq/2G0-  2a/(al-4)[Uln(r/c)]-(a~3)/(o-2)a(R/r)(a"2)  (3.4.18b) 

-(a+1)/(a+2)l(R/r)(ct+2)  , 

a  /2G„-  1/(a-2)>2a/(a2-4)ln(r/c)-(a-3)/(a-2)l(R/r)(a-2)  (3.4. 18c) 

ww 

+1/(a+2)2(R/r)(a>2)  , 

where  the  constant  c  will  be  determined  from  the  boundary  conditions. 

Suppose  that  the  undeformed  surfaces  of  the  block  a,-0,  a2-H,  b-L/2, 

and  c-1/2  in  Fig. 2  map  into  the  deformed  surfaces  r-r,t  r-r2,  and  3-90.  The 

following  geometrical  restrictions  must  hold: 

(a)  R-(ri-r2,)/(aO; 

(b)  90-L/2R  and  in  the  case  in  which 

the  beam  is  bent  as  a  circle,  R/H-L/(2Hir) .  (3.4.19) 
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T  -  2a/(a*-J*)C(n2/2)m(n*/d)-(n?/2)ln(nl/d)  +  (ni-nf)/'0 

>(a-3)/(ct-2)2/(a-M)[(n/n3)(a''2)ni-  (n/ n, )  (a"2)  n?  ]  (3.4.21) 

♦(a+1 )/(a*2)2/a  C(n/na) ^a*2^n|-  ( n/ n, )  ^a+2^ n?  ] 

S  -  n*gj-  nig, 

n  -  (n|  -  n!  )/2 
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ni-  r>/ H  ,  n,«  r*/H  ,  n  -  R/H  ,  T  -  M/(2G0H*),  d  -  c/H 

S  -  F/(2G0H)  ,  g,-  fl/(2G»),  and  g*-  f,/(2G0) 


(3 • 4.22) 


The  problem  of  Interest  here  is  the  one  in  which  the  radial  boundary 
: stresses  f,  and  f2  are  equal  to  zero.  The  solution  then  reduces  to  assigning 

a  value  to  R,  the  radius  of  curvature  of  the  neutral  axis.  Eq.  *s(3. 4.21 ) 
then  determines  the  deformed  radii  r,  and  r2,  the  constant  of  integration  c, 

the  moment  M  and  the  end  force  F  which  is  zero  in  this  case. 


Specific  cases 


case  (i):  A  very  thin  beam. 

As  the  first  example,  consider  a  pure  bending  of  very  thin  beam.  Let  £ 
be  the  distance  from  the  neutral  axis.  Th-n  i--R+£  ,  and  £/R<<l  .  Employing 
the  Taylor  expansion, 

(R/r)(at"2-  1-(a-2)U/R) ;  (R/r)(a+2-  Ma+2)U/R) 


(3-4.23) 


The  non-zero  stresses  (3.4.15)  become 


<7rr-  -P-2G,£/R  ;  ff00-  -P*2G0£/R  ;  -P 


(3.4.24) 

Satisfaction  of  the  equilibrium  eq. (3.4.16)  gives  the  pressure  as 


P  -  -2G0£/R  *  c 


(3.4.25) 


introducing  (3.4.25)  into  (3.4.24), 


a  -  -c  ,  <j  -  4G05/R  -c  ,  o  -  2G05/R  -c  .  (3.4.26) 

rr  99  ww 

The  boundary  condition,  orr“°  .  at  £-£i,£2  gives  the  integration  constant 

c-0;  and  the  condition  of  zero  resultant  force  at  the  ends  of  8-dlrection 
yields  £i*-'£2-h/2  .where  h  is  the  deformed  thickness  of  the  beam.  Then, 

o  -  4G0£/R  ,  a  -  2Ga&/R  ,  M  -G„h,/(3R).  (3.4.27) 

o  “  WW 

The  results  here  are  the  same  as  those  in  the  elastic  deformation  of  pure 
bending  except  the  stress  in  w-direction  which  is  required  to  produce  plane 
strain  conditions  in  the  presence  of  incompressibility.  Evidently  the 
strains  in  this  case  are  extremely  small  resulting  in  an  elastic  solution. 

Case(li):  A  thick  beam. 

Now  consider  an  example  of  pure  bending  of  a  thick  beam.  In  thi»  case 
the  zero  resultant  force  F-0  and  f,-f2-0.  Therefore,  gi-g2-0.  thus, 

2a/(ai-4)ln(nl/d)*1/(a-2)2(n/ni)(a_2)  +  1/(a*2)2(n/n1)(<i+2)-  0 

2a/(a2-4)ln(n2/d)  +  1/(a-2)2(n/ni)  (a~2)-M/(a+2)2(n/n2) (a+2)-  0 

n  -  (n|  -  n?  )/2  (3.4.28) 

T  -  2a/(a2-4)[(n2/2)ln(n2/d)-(n?/2)ln(ni/d)+(n2-n2)/4] 

^(a-3)/(a-2)2/(a-4)C(n/n2)(a‘2)n!-  (n/ni)(a”2)n?  ] 

+  (a+1 )/(a+2) 2/a  [(n/n2)^a  2^ni_  (n/ni)^a  2^n?  ] 


Eq. (3*4.28)  is  a  set  of  non-linear  equations.  The  following  procedure  is 
used  to  get  some  numerical  results  for  this  set  of  equations. 

(1) .  Eliminate  the  constant  d  from  the  first  two  equations  of  (3. 4.28)  to 

get  one  equation  involving  unknown  n, ,  n*.  and  r\. 

(2) .  Use  the  Newton  iteration  to  solve  Hi  and  n2  for  given  n  from  the 

resulting  equation  from  step  (1)  and  the  third  equation  in  (3.4.28), 
which  are, 

2a/(o2-4)ln(n,)>1/(a-2)2(n/hi)(a“2)>1/(a^2)i(n/rh)(a+2)- 

2a/(a2-4)ln(n1)>1/(a-2)2(n/n2)(a"2)  +  l/(a+2)2(n/r12)(a+2) 

n  -  (n|  -  n?  )/2  (3.4.29) 

(3) .  Obtain  the  constant  d  for  corresponding  n,  hi,  and  r\z  from  either  one 

of  the  first  two  equations  in  (3.4.28). 

4).  Obtain  T  from  the  last  equation  in  (3.4.28)  for  corresponding  n ,  n i  » 

Hj ,  and  d . 

Before  formulating  the  numerical  solution  scheme,  the  existence  and  the 
uniqueness  of  eq.(3.4.29)  is  first  discussed.  Look  at  the  following 
function, 

f(x)-2a/(a2-4)lnx+1/(a-2)2(n/x)(a“2)+1/(a+2)2(n/x)(o+2)  (3.4.30) 

where  x  may  be  n»  and  n».  For  given  n,  when  x  approaches  zero,  f(x)  goes  to 

Infinity.  When  x  tends  to  infinity,  f(x)  goes  to  infinity.  f(x)  has  one  and 

only  one  minimum  at  x  approximately  equal  to  n,  since  a  is  a  large  number 

m 

ranging  from  100  to  200,  x  -n  is  the  root  of  following  equation, 


r(x)'-C2o/(a2-4)-(n/xm)(a'2)/(a-2)  +  (n/X[n)(a+2)/(a+2)]/xm-0  (3.4.31) 

Fig.  15  depicts  the  character  of  function  f(x).  It  can  be  seen  that  for  a 
value  of  f(x)t  there  are  two  roots( except  the  case  in  which  the  minimum  is 
also  the  root)  corresponding  to  nt  and  n2.  One-one  relation  I  is  shown  in 

Fig. 16  from  the  first  equation  (3-4.29).  On  the  other  hand  another  one-one 
relation  II  is  also  shown  in  Fig.  16  for  positive  hi  and  n2  from  the  second 

_ _  equation  of  (3-4.29).  The  intersect  of  I  and  II  gives  the  solution  m  and  n2 

for  given  n. 

The  Newton  iteration  scheme  [56]  for  (3-4.29)  can  be  written  as 


hj 

hi 

30j/3ni 

3$i/3n2 

-1 

n2 

i  +  l 

hi 

i 

3$a/3n i 

k 

3$2/3n2 

i 

♦  * 

(3.4.32) 


i  , 


where  and  <j>2  are 

-  2a/(a2-4)in(hi/nl)  +  [(n/nl)  2^-(n/n2)^a  2^]/(a-2)*  + 

♦  C(n/n1)(<S>2)-(n/ni)(a+2)]/(o*2)3  (3.4.33) 

and  -  2n  -  n|  +  n? 

Numerical  results  were  obtained  after  setting  a-200.  Fig. 17  shows  the 
variation  of  moment  with  the  curvature.  There  is  a  softening  effect  in  the 
sense  that  the  moment  decreases  with  increase  in  curvature  after  a  critical 
value  of  the  curvature  has  been  reached.  However  this  is  a  geometric 
Phenomenon  since  all  stresses  Increase  raonotonically  with  curvature.  Fig. 18 
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shows  the  variation  of  o 


with  curvature,  at  the  neutral  axis  while  in 


pig. 19  and  20  we  show  aQ9  at  r-rj  and  r-r2  .  A3  the  curvature  increases  the 
neutral  axis  moves  closer  to  the  compressive  surface  of  the  beam  thereby 

•  £ 

requiring  a  large  compressive  tangential  stress  to  maintain  moment 

fi? 

r  equilibrium.  The  material  can  sustain  a  large  compressive  stress  by  virtue 

of  its  incompressibility.  Fig. a  21-23  show  the  distribution  of  c  over  the 

0  0 

!  Jrv 

cross-section  of  the  the  beam  for  values  of  R/H  of  50,  5,  and  1.  The  stress 

a  is  zero  at  the  neutral  axis  as  expected.  Fig.s  24-26  show  the 
99 

distribution  of  a  over  the  depth  of  the  beam  for  values  of  R/H  of  50,  5, 
rr 

and  1  .  Observe  that  o  is  about  two  orders  of  magnitude  smaller  than  aQQ 
and  reaches  maximum  value  at  the  neutral  axi3. 


3.5.  Symmetrical  expansion  of  a  thick  spherical  shell 

Identify  the  curvilinear  coordinate  systems  in  undeformed  body  and 
in  deformed  body  with  the  system  of  the  spherical  coordinates  (R,$,Q)  and 

(r,$,9),  respectively.  The  deformation  of  the  symmetrical  expansion  of  a 
thick  spherical  shell  can  be  described  as 


r  -  r(R,t)  ,  <*•'$,  9  •  9  . 


(3.5.D 


In  the  material  system  metric  tensors  is 


1  0  0 


|C  .|-R"3Ui!«  (3.5.2) 

BP 


0  0  (Rsin$)2 
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and  in  spatial  system 


[gag]  -Or*  0 


[g^l-r'sin1*  .  (3.5.3: 


0  0  (rsin$)2 


-«*.■  The  deformation  gradient  matrix  is  given  below, 
’r.H-  \  r'  0  0  1 


CF]  -  0  1  0 


(3-2.4) 


[0  0  1  |  , 

where  r'-dr/dR.  The  Cauchy-Green  deformation  tensor  i3 
I  r'  0  0  1 


[C]  -  0  r2  0 


C  -  r'  2r*3in2(j> 


(3.5.5) 


[  0  0  (rsln$)2J 

For  incompressible  material, 


r^r'-R"  ,  i.e., 


r2dr-±R2dR 


(3-5.6) 


Positive  sign  in  eq.(3.5.6)  is  taken  for  expansion  of  the  sphere.  After  an 
integration,  r1-  R*+A(t).  Therefore,  the  deformation  (3-5.1)  becomes 


r  -  [R2+A(t)]1/3, 


♦  *  ♦  * 


9-0 


(3.5.7) 


Consequently,  the  Cauchy-Green  deformation  tensor  becomes 


RVr*  0 


[C]  -  0  r2  0 


C  -  r '  2r‘,sin2ii 


(3.5.8) 


0  (rsin$)2 


-  50  - 


then  the  inverse  of  the  Cauchy-Green  tensor  is 
rVR*  0  0 

[C]  '  -  I  0  1/r*  o 


-1 


0  0  1/(rsin$)2 

while  the  increment  is 


d[C]  - 


-4(R*/rs)dr  0 


2rdr 


(3.5.9) 


(3.5.10) 


0  0  2rsin2$drj 

The  intrinsic  time  can  be  obtained  by  substitution  of  eq.(3.5.9), 
eq. (3.5.10)  into  eq. (2.2.22)  as 

(dc)2-(d r/r)2  ,  (3.5.11) 

where  p2  has  been  set  to  1/6.  For  monotonic  expansion,  d^-dr/r.  After  an 


integration,  ;-ln(kr),  where  k  is  a  constant  of  integration.  Since  c«0,  for 
r-R,  hence  kR-1 ,  i.e.,  k-l/R.  Neglecting  the  hardening  and  softening  effect, 
dz  and  z  are  obtained  as  follows, 

dz-dr/r  ,  z-ln(r/R)  .  (3-5.12) 

With  the  h  .*lp  of  eq.  (3.5.12)  d[C]/dz  is  given  below, 


d[C]/dz  - 


-4(RVr*)  0 


2rs 


0  2r2sin2$ 


(3.5.13) 


A 
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Taking  G-G0e~aZ,  the  following  stresses  are  derived  by  use  of  the 


constitutive  equation  (2.2.25). 

4G«  o-4 

V  -  ur)n-  £’-<»"•>  3) 

2G0 


rf 

>  *  .'tf*  - 

trr' 

- 

;  s&j 

t'.-iqjqw 

:  * 

X 

T00 

^r. 

■»--  and 

i  4-' 

TR0 

(a*  2) 


[1-(R/r)a+2]}  , 


2G# 

(a+2) 


Cl-(R/r)a+2]} 


(3.5.14) 


i;  The  components  of  stresses  in  spatial  curvilinear  system  are 

4G«  _ n 


T  -  -P  " 
rr 


T^T 


rr  [l-Ol/r)*""  ]  . 


V  rll-*p  *  T^5T  , 


700-  r2sin2*  {-P  ♦ 


[3J5T  Ci-(R/D«  2]| 


(3-5.15) 


a«d  T  -  T.a-  T  -  0 
re  *e  <pr 

The  physical  components  of  the  stresses  are 


°rr"  "P  "  (a-4)  Cl-(R/r)  3 


o  a+  p 

o  -  aoa-  ~P  ♦  t — rr  [1-(R/r)a  *] 
$4  99  (a+2) 


(3-5.16) 


“d  ar9"  %e"  V“  ° 

For  the  symmetrical  expansion,  the  equilibrium  equations  in  spherical  system 
with  the  absence  of  the  body  forces  are 
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tfpr • 


where  R,  and  r,  are  inner  .  adii  of  undeformed  and  deformed  sphere, 

respectively.  With  the  help  of  (3*5.19),  the  following  non-vanishing 
stresses  are  obtained, 

a  -  4G,Q(R/r)  -  o 
rr 

j|Q  2q 

%♦*  »rr*  CMa/D'-V-j^y  (3.5.20) 

Consider  the  following  boundary  conditions  in  the  resent  problem, 

(1)  «rr"pi  »  at  r-r,  . 

(ii)  o  -P2  ,  at  r-r,  .  (3.5.21) 

rr 
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;3 

,4^v  ^ 


Upon  satisfying  the  condition  (i)  in  eq.(3.5.2l),  c— Pt  .  Applying  the 


condition  (ii), 


Pi  -  Pi  -  4GaQ(R ,/r,) 


(3.5.22) 


From  the  geometrical  condition  r*-R**a(t), 
rj-  R!  -  r*  -  R>  -  r*  -  R‘  , 


(3.5.23) 


ss 


where  Ra  and  r2  are  the  outer  radii  of  the  undeformed  and  deformed  sphere 

respectively.  After  a  few  steps  of  simple  algebra  manipulation,  the 
^  following  relation  is  derived, 


CK*-1  ♦  l/tRj/rj)*]173 


(3.5.24) 


where  K-R/R, .  Eq.s  (3.5.22),  (3.5.24),  (3. 5. 20)  and  (3.5.19a)  with  the  known 

P»*  P*»  P 1  and  Ra  give  the  complete  solution.  It  can  be  seen  from  (3.5.22) 

that  when  P,-PJt  Q(Ra/ra)»0.  Since  the  integrand  of  eq . (3. 5.19a)  is 

monotonic  increasing  function  it  follows  that  R1/r1-Rl/ra.  This  is  true  if 

only  if  Rj/r^  Ra/ra-l.  This  is  the  condition  under  which  no  deformation 

takes  place.  This  agrees  with  the  Incompressibility  assumption.  Some 
numerical  solutions  with  R2/Ri*»1.5  and  a-200  are  given  in  Fig. 28-30.  Fig. 28 

shows  the  variation  of  (Pa-  P,)  with  respect  to  rj/R,.  The  former  decreases 

after  a  critical  value  of  r./R,.  Fig. 29  shows  the  distribution  of  the  0 

rr 

along  r-direction  for  a  case  P,/ 4G9 -0.0051  and  Pa-0  while  Fig.  30  shows  the 
distribution  of  cQ0  and  a  along  r-directlon  for  the  same  case. 
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In  this  chapter,  the  closed  form  solutions  for  five  special  problems 
were  obtained,  in  which  the  prescribed  deformation  field  has  one 
undetermined  parameter.  But  for  more  complicated  problems  it  can  be  expected 
that  the  numerical  solutions  are  necessary.  In  next  chapter,  the  numerical 
method  with  finite  element  technique  for  the  boundary  value  problem  win  be 
formulated  for  both  incompressible  and  compressible  plasticity. 


Fig- 6  The  deformation  profile 
for  various  parameter  5 


FI 5. 7  The  relation  between 

aeformatlon  u,  and  the 

tne  parameter  S  (r-0). 


Fig. 3  The  distribution  of  c 

and  a  over  the  cross 
ww 


Fig. 9  Torsion  of  a  bar 


Fig.  10  The  relation  between 
torque  and  twist. 
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1 


f(x)-2a/(al-M)lnx-1/(a-2)*(n/x)(a"2)^1/((a+2)l(n/x) 
Fig.  15  Typical  plot  of  fCx) 
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(a*  2) 


I- 2a/(aJ-«)ln(n,/na)*l/(a-2)1 II-[(n/n»)(a“2  >-(n/na)(a"2  ]1 

♦  l/(a-2)*C(a/T,l)(a*'2  )-(n/iia)(o'>2  }] 

II- (ni-n?)/2 

Fig. 16  The  plot  of  function  I  and  II 
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Fig. 18  Stress  9 _  at  neutral 

rr 

axis  vs.  tne  curvature 
for  a«200. 


Fig. 20  Tha  variation  of  Stress  ag9 
at  r-r,  for  a- 200 
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MJ  **J  9U  3U  S'.J 

r/H 

Fig. 21  The  distribution  of  o80 

over  cross  section  for 
R/H -50  and  a- 200 


■ft**  ^ ^ 

4  i  i 

r/H 

Fig. 22  The  distribution  of  efte 

over  cross  section  for 
R/H -5  and  a- 200 
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Fig. 2*  The  distribution  of  b 

r . 

over  cross  section  for 
R/H-50  and  a-200 


CHAPTER  4  -  ENDOCHHONIC  CONSTITUTIVE  RELATION  OF 


X 


COMPRESSIBLE  PLASTICITY  WITH  FINITE  DEFORMATION 

The  assumption  of  plastic  incompressibility,  which  is  pivotal  in  the 
development  of  the  classical  theory  of  plasticity,  is  an  approximation  and 
simplification  to  the  real  situation.  For  some  rubber-like  materials  this 
assumption  may  be  a  good  approximation,  but  it  is  not  a  universal  rule,  in 
reference  [57]  the  assumption  was  evaluated  through  a  series  of  3imple 
tension  tests  on  some  very  important  materials  such  as  aluminum,  copper  and 
low  carbon  steel.  The  experimental  results  showed  that  those  materials  are 
plastically  compressible  and  that  the  compressibility  Increases  with 
straining  in  the  plastic  region.  Hence  developing  a  theory  for  compressible 
plasticity  is  very  Important  task.  In  this  chapter,  the  plastic 
complressibility  will  be  introduced  into  the  endochronic  constitutive 
equation.  To  do  this  the  free  energy  density  in  the  thermodynamic 
formulation  will  be  first  modified.  Then  the  method  used  for  developing  the 
constitutive  relation  for  incompressible  plasticity  in  the  Chapter  2  will  be 
used  to  get  the  constitutive  equation  in  the  presence  of  compressible 
plastic  deformation. 

**.1  Development  of  the  Constitutive  Equation 

In  chapter  2  a  constant  term  in  the  free  energy  form  was  assumed.  It 

is  noticed  that  this  assumption  led  to  the  constitutive  equation  of 
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f 


incompressible  plasticity.  Now  the  constant  term  is  replaced  by  a  function 
of  It  >  the  third  invariant  of  the  Cauchy-Green  deformation  tensor.  Hence, 


the  free  energy  density  becomes 

*  ,»  it  \  „(**)  (r)  ..-..(r)  ( r )  (r) 

f  -  fa(I,)  +  A  qli  ♦  1/2C  q^j  qiJ 


(4.1.1) 


Along  the  line  of  the  development  of  the  constitutive  relation  for  the 
incompressible  material,  it  can  be  seen  that  with  this  modification  to  the 
free  energy  density  the  covariant  components  of  q^  in  the  material  frame 

remain  of  the  form 


1(P). 

*afl 


(r) 


(r) 

n  _  A  ,z  -a  (z-z')2  ,  , ,  .  , 

T7T;0e  p  C„B(I  > 


,(P) 


(4.1.2) 


where  a 


r  b(r] 


Fo"  the  constitutive  relation  now  the  contravarlant  components  of  the 
stress,  upon  use  of  eq.(2.l.8)  in  conjunction  with  eq. (4.1.1),  are 


_aB  o  r?T  _ 1  ra0  .  r(r)  (r)a  (r)B’  n 

7  *  p,C2I>  3i7  c  2(a  q(r)  C  q  e*q  “  )] 


B 


'(r) 

or  the  covariant  components  of  the  stress 

r  .  £  r?T  —  r  ?u(r)n(rItr(r)n(r)n(r)8'n 
Tafl  1C21*  317  CaB  2(A  qaB  *  qaS-q  B  )] 


(4. 1.3) 


(4.1.4) 


Substitution  q.(4.1.2)  into  eq.(4.i,4)  gives  the  following  relation 
for  the  covariant  components  of  the  stress 
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T  -  £  [21,  r—  C  +  /*G( z-z'  )C  Cz*  )  dz' 
a8  pa  31,  a8  0  aB 


dC  dCa 

-Cuv^;^G(2z-z1-zl)-3iSi  _±i  dz,dz2]  . 


(4.1.5) 


Once  again  ignore  the  small  quantity  of  the  double  integral  in  eq.(4.1.5) 
and  get  the  constitutive  equation  as 

3tt°«8 


In  the  matrix  form 


[T]  -  J  {21,  gf7  CC3  *  /qG(z-z’  )[C(z'  )]  dz'  } 


(4.1 ,6a) 


.  (r) 

G(z)  -2  l  e~arZ 

r  CT  J 


(4.1.6b) 


When  the  condition  of  incompressibility  is  imposed,  |c  .1-1,  the  term 

a8 

3f,/3I,  is  indefinite.  If  the  term  is  set  to  be  -P  (the  factor  21,  is  also 

absorbed),  eq.(2.1.25)  is  again  derived,  which  is  the  constitutive  relation 
for  incompressible  plasticity. 


T  o"  “p  c  a  +  /?G( z-z ' )C  Cz')  dz’ 
aB  aB  0  aB 


(2.1.25) 


For  Intrinsic  time  scale,  d^  is  no  longer  vanishing  under  a 
compressible  plastic  deformation,  thus 


(  3t  >*  -  '’"Vjj  *  >*Vij 

Eq. (2.1.15),  after  using  eq .(2.2.16),  becomes 


(2.1.15) 
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<  if  >•  -  P,/»  «=;;  csa)*.  p,  /«<=;!  £as  e;;  cnY 


(4.1 .7) 


<C  “.a1**  P!  /4  °ri  dc  oa  csn  dc„Y 


(4.1 .8) 


In  matrix  form 


dc2  -  Pj/4  {tr[B]}2+  p2/4tr[D]  , 


(4.1.9) 


where 


[D]  -  [cf’dCClCcf’dCC] 


(4.1.9a) 


[B]  -  CC]  dCC] 


(4.1.9b) 


An  intrinsic  time  scale  z  is  still  defined  to  account  for  the  hardening 
or  softening  character  of  a  material, 


rZ  dt’ 

f ( c' 


(4.1.10) 


where  f(;)  is  a  positive  function  of  c.  If  no  hardening  (softening)  takes 
place  f( 0-1  • 

In  the  curvilinear  system  as  described  in  section  2  of  Chapter  2,  the 
curvilinear  components  of  the  Cauchy-Green  deformation  tensor  were  given  by 


39k  30 

C  -  g  — -  — - 

uv  kl  3*“  3*v 


(2.2.8) 


and  the  curvilinear  components  of  stress  T  .  by 

ap 


T  T  3»a  3»8 
uv  ’  08  3XU  3XV 


(2.2.23) 
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Upon  use  of  eq. (2.2.8)  an  alternative  form  of  the  intrinsic  time  can  be 
derived, 

(  )»  -  »./*«;;  «,„>**  P.  /*  si  «ae  «nT 

or,  in  matrix  form, 

(  dc  )a  -  Pj/4  ltr[B])*>p2/4  tr[D]  (4.1.12) 


where 


ft]  -  [C]_1d[C]  CC]_1d[C] 


(4.1.12a) 


and 


[B]  -  [cf’dCC]  (4.1.12b) 

The  constitutive  eq. (4.1.6)#  by  use  of  eq. (2.2.23)#  becomes 


3f  1 

x  -  [ 21 ,  rr~  C  ♦  /?G(*-z')C  (z')  dz'] 

uv  Po  *  31,  uv  0  yv 


(4.1.13) 


or 

gy  “ 

[t]  -  £  (21.  xrA  EC]  ♦  /^G(z-z’)[C(z')]  dz*}  .  (4.1.13a) 

Po  31| 

4.2  Discussion  of  the  Function  ?<, 

In  last  section  the  function  %  appears  as  a  functional  form  in  the 

constitutive  equation.  How  to  choose  the  function  so  that  it  can  describe 
the  compressible  behavior  of  a  material  well  needs  much  effort.  In  this 
paper  a  particular  form  of  given  by  eq.  (4.2.1)  is  taken  to  demonstrate 

the  application  of  the  constitutive  equation: 


f0  -  U/8)[ln(I,)]1  ,  (4.2.1) 
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where  X  is  a  material  constant. 


It  is  now  shown  that  the  constitutive  equation  with  the  assumed  f0  is 

consistent  with  elasticity  theory  under  infinitesimal  deformation. 
Substituting  eq. (4.2.1)  into  the  constitutive  equation  (4.1.6), 

T  -  p/p0{Xln(J)C  ♦  G(z-z’)C  Cz')dz*}  (4.2.2) 

op  up  up 

where  J  is  the  Jacobian  defined  as 

J-l!/2-| 3Y./3Xa|  .  (4.2.3) 

When  a  deformation  is  very  small,  p/p0-1,  C  -6  and  without  lose  any 

up  up 

generality 

J-X,X2X,-(1+en)(1+e22)(i  «■£,,) 

*1  +  (ex  j  *e2  2 ♦e , 2 )*0( e* )  »  1+e  ,  (4.2.4) 

where  X ^  and  ^  (i  — 1 ,2,3)  are  principal  values  of  the  deformation  gradient 

matrix  and  the  strain  tensor,  respectively.  Hence, 

ln(J)  -  ln(1+e^)  •  e^+0(e^)  ■  •  (4.2.5) 

Eq. (4.2.5)  represents  the  dilation  of  a  material  for  Infinitesimal 
deformation.  With  the  help  of  eq. (4.2.5),  the  Hooke's  law  is  got, 

TuS-  Xeli5uS+  2G(0)eaS  (J4-2-6) 

where  X  is  a  Lame  constant  and  G(Q)  the  shear  modulus. 

In  next  chapter  the  constitutive  relation  given  in  eq. (4.2.2)  will  be 
used  to  analyze  the  problems  by  means  of  the  finite  element  technique. 
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CHAPTER  5  -  NUMERICAL  SOLUTION 


A  numerical  algorithm  baaed  on  the  finite  element  method  of  analysis 
of  the  boundary  value  problem  is  developed  in  this  chapter.  The  formulation 
of  the  finite  element  equations  is  referred  to  the  material  system,  namely, 
the  Lagrangian  formulation.  The  finite  element  formulations  for  both 
compressible  and  incompressible  material  are  developed.  A  computer  code  is 
established  to  solve  a  plane  strain  boundary  value  problem  under  a 
compressible- plastic  deformation  by  use  of  the  linear  triangular  element. 
Finally,  the  solutions  of  upsetting  of  a  block  (a  forging  process)  are 
presented. 


5.1  Fundamental  Equation  (the  Principle  of  Virtual  Work) 

Equations  of  equilibrium,  when  body  forces  are  absent,  in  rectangular 
coordinates  of  a  material  system  can  be  shown  to  be 

(JT®  )  -  0  (5.1.1) 

1  .a 

where  (JT01^)  is  the  stress  per  unit  undeformed  area.  The  mixed  stress  T0^  , 


which  is  the  projection  of  the  force  T  along  i-direction  in  spatial  system, 
is  defined  as 


T  i"  °ij  Xj 


9lj"  T  lyj,a 


(5.1.2) 


where  x®-  3xa/3yj ,  is  the  Cauchy  stress,  (  )f(j  represents  3()/3xa. 
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After  defining  U  to  be  the  displacement  and  u^  the  components  of  the 

displacement  and  5u^  the  components  of  the  virtual  displacement, 

Multiplication  of  Su^^  to  the  equilibrium  equation  (5.1.1)  and  then 

integration  of  the  resultant  equation  over  the  undeformed  volume  give  the 
following 

/„  (JTa. ),  $u  dV4  -  0  .  (5.1.3) 

The  application  of  Green's  theorem  to  (5.1.3)  yields  the  equation  of  the 
principle  of  virtual  work  under  large  deformation  as 

/„  (JT°. )5u.  dV„  -  /  ( JTa. )u  dS„  ,  (5.1.4) 

Vq  X  1  9  Cl  WQ  11 

where  the  surface  integration  of  the  right  side  of  (5.1.4)  i3  over  the 
undeformed  surface.  Equation  (5.1.4)  will  be  the  basis  of  the  finite  element 
formulation. 

5.2  Incremental  Form  of  the  Endochronic  Constitutive  Equation 

In  Chapter  4  the  endochronic  constitutive  equation  for  compressible 
material  was  derived  as 

JTag-  Uog(J)Cag+  /* G(z-z')C  jjZ')  dz'  .  (4.2.4) 

The  incremental  form,  with  respect  to  the  intrinsic  time  scale,  can  be 

obtained  from  the  intrinsic  time  derivative  of  (4.2.4)  as 

*  _  * 

d(JT  )/dz  -  (X/J)JC  ♦  G  C  ♦  dQ  „ 
aS  aS  a6  a6 


(5.2.1) 


or 


d(JT  J  -  (X/J)  dJ  C  ♦  G  dC  .  «■  dQ  dz 

aB  aB  aB  a0 


(5-2.2) 


where 


and 


G  -  Xlog(J)+G(0) 


X  dG(z-z’)  dCot8 
dQa6*  ^  - - dF 

For  an  Incompressible  material. 


(5.2.3) 


“(JTa8>  '  -  UP)C<,B*  0  ‘“.s'  *  “.s“ 


(5.2.4) 


where  G  -  G(0)-P  and  dQ  .  remains  the  same  as  in  (5.2.3). 

ao 


5.3  Formulation  of  the  Finite  Element  Equations 

The  finite  element  approximation  is  developed  from  a  displacement 
assumption  within  each  element,  which  gives  the  displacements  at  any  point 
within  the  element  aw  a  linear  combination  of  the  displacements  at  a  finite 
number  of  nodes  of  the  element,  the  coefficients  being  constants  or 
functions  of  the  position  within  the  element.  In  this  section  we  adopt  such 
a  displacement  assumption  in  general  form  using  material  coordinates,  i.e., 
a  Lagrangian  formulation. 

Assuming  the  displacement  within  any  element  in  the  form 

ut-  ♦®(x°.  xa)qm  ,  (5.3.1) 
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where  q  denotes  the  nodal  displacements;  $  is  the  transformation  function 
m 

which  gives  the  displacements  at  any  point  within  the  element;  xa  are  the 
nodal  coordinates. 

The  following  equations  are  some  derived  geometric  relations  which  will 
be  used  in  the  finite  element  formulation. 
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(i»)  Since  we  have  the  relation 


yK,a  XJ  " 


(5.3.8) 


x?-  [y4  J' 


(5.3.9) 


(5)  A  differentiation  of  both  sides  of  eq. (5.3.8)  gives 


!l(yk,»,x?  *  yk,a  “<XJ)  *  0  ' 


(5.3.10) 


Then  what  follows  is 

4xk  ■  •  x“  4  ari.r  -  x“  4  ♦Tb'V  ■  xi  •  (5-3-"> 

where 


XK  ♦ifl  “  (defined) 


(5.3.12) 


(6)  The  Cauchy-Green  deformation  tensor  and  its  incremental  form  are 

Cat~  yk,<l  yK,  8 


dCaS-  (Cyk,a'^SyK,aMV  ^Okj’Wl .  •  (5-3',3) 


where 


♦(ik>-  1/2(*ik+  ♦ki)  (deflned)  • 


(5.3* 1  *•) 


(7)  The  Jacobian  and  its  incremental  form  are 

J  -  |y.  .1 
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dJ  -  Jx“  *“adqm- 


(5-3-15) 


Having  the  assumed  displacement,  we  get  the  equilibrium  equation  for  an 
element  by  direct  substitution  of  (5-3-1)  and  (5-3-2)  into  the  equation  of 
the  principle  of  virtual  work  (5.1.4) 


V^Xa5^  VJT>5VadS- 


(5-3-16) 


Since  the  virtual  displacement  Sq^  (m-1,2 . )  are  Independent  of  each 

other,  the  following  m  equations  are  obtained 


VJT°Xadv*-  V 


(5-3-17) 


where 


fm  '  VJTV<VS- 


(5.3.18) 


The  Incremental  form  of  (5-3-17)  can  be  written  as 


V<JTV  alv,x 


(5.3.19) 


wnere 


d(f)  „  -  /<.  d(JT°.)  dS, 
n  ex  So  i  la 


(5-3-20) 


From  the  tensor  transformation 


„  .  a  J  T 
°ki  k  i  Tafl 


(5-3-21) 
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and 


Finally,  (5.3.19)  in  conjunction  with  (5.3.27)  results  In  the  following,  a 
set  of  j  linear  equations  In  dq^  , 


;v.c**Mk*  25,;ik)*u'2,!k3i'*ik!J03i)'*"i'»L!Jok3nav" 


it*LdQu<aa',‘ 


(5.3.28) 


Alternatively,  (5.3.28)  can  be  written  in  matrix  form  as 

[k.  ]  dq  -  d(r.)  *  d(f.)  -  d(f(.t))  , 

Jm  m  J  ex  J  p  j 

where  the  element  stiffness  matrix  k.  is 

Jm 

V  V'Okk*  26*Uk),‘,ik‘2*,(l<a)tJil<(J<’sl,'*"i'*Jik(J'’kS,:illv" 


(5.3.29) 


and  the  incremental  pseudo-force  vector 

a<Vp  ‘  '  Vu'Wki‘te'JV"  • 


(5.3.30) 


where  the  integration  is  performed  over  the  undeformed  area  of  an  element. 

Considering  the  equilibrium  of  the  entire  structure,  we  obtain  the 
structure  stiffness  matrix  in  the  form 

[K]  {dq}  -  {dF}  (5.3.3D 

where  [K]  is  global  stiffness  matrix,  {dq}  the  total  incremental  nodal 
displacements,  {dF}  the  total  Incremental  nodal  force  which  consists  of 
applied  force  and  the  incremental  plastic  pseydo-force  calculated  using  the 
second  equation  of  (5.3-30). 

To  complete  this  section,  the  formulation  for  incompressible  materials 
is  derived  through  the  mixed  finite  element  method  as  follows.  The 
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displacement  and  the  pressure  are  taken  as  the  primary  variables.  The 
displacement  within  any  element  remains  in  the  form  of  (5.3.1), 

ui”  *i(V  -°)qm  (5.3.D 

and  the  pressure  within  the  element  is  approximated  as 

P  -  pn(xa,xa)pn  ,  (5.3.32) 

where  p^  is  the  nodal  pressure.  The  degree  of  the  interpolation  function  pn 

for  the  pressure  should  be  one  order  less  than  the  degree  of  the  function  41 
for  the  displacements  in  order  to  have  consistency  of  approximation  for  the 
displacements  and  the  pressure.  This  will  be  seen  clearly  from  the  following 
derivation  because  P  is  one  order  less  differentiable  than  u^ .  The 

incremental  form  of  P  is 

dP  -  pn(xot,x°)dpn  .  (5.3.33) 

When  the  incremental  form  of  constitutive  equation  for  incompressible 
material  (5.2.4)  is  used  in  the  last  term  of  eq. (5.3.24)  with  J-1  (hereafter 
in  this  section) , 


A  \  A  d<TY8>  ■  A  A  xf  [-')P  CYS*  5  a(CY6>*  '“>•,.<*3 


(5.3.34) 


89  - 


Substituting  (5.3-34)  into  (5-3.24), 


d(T°1)-x“(C2G^ll<r^ks)(a3i)-^i(0i<s)]dqm-  dQ^dz}  .  (5-3.35) 


Finally,  (5-3- ^9)  in  conjunction  with  (5-3.35)  results  in  the  following 

/v.C25*Uk)*3lK'2*(k3)*Jlk(J“sl)'*"l*lK<j0l<S):llIV"d'1.'  V0"»nau‘apn 


alfJ>ex-  ;v„*uaQkl“av-  ■  (5-3-3«> 


In  matrix  form  eq. (5-3-26)  becomes 


*  [‘■jn^n  -  dlfJ>ex  *  d(fjV  d“j  >  ' 


(5.3.37) 


where  the  element  stiffness  matrix  k.  Is 

Jn 


V  ;v„[25*“lk)*Jik'2*(k3)*JlkUl,31>‘*si*Jlk(J,’k3):ldV-  ■ 


L  -  -  /  pV  dV„ 
jn  JV0P  viL  0 


(5.3-38) 


d(Vp  “  -  /V0^kdQkidZdV* 


Now  we  have  m  equations  but  m+n  unknowns  (m  for  displacement  and  n  for 

pressure).  The  constraints  of  incompressibility,  i.e.,  J-1  or  j-0  must 
considered.  From  equation  (5-3-15), 


.  dq  -  0 
11  m 


(5-3-39) 
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If  we  multiply  p  to  both  aides  of  eq.(5.3.39)  and  then  integrate  the 
resultant  equation  over  the  undeformed  element,  we  can  obtain  n  additional 
equations  given  by  eq.(5.3.40) 


“"a  ‘ 0 


(5.3.40) 


In  the  matrix  form,  (5.3*40)  becomes 

[L  ]dq  -0 
nm  m 

The  simultaneous  equations  of  q's  and  p's  can  be  written  as 


(5.3.41) 


k , 
jm 

S)n 

dqm 

L 

nm 

0 

dPn 

k 

m 

0 

(5.3.42) 


For  an  entire  structure  we  have  the  equation  in  the  following  form. 


*NN  S(M 
Sin 


'dqM 

Jt) 

dF 

rN 

dPN 

» 

0 

(5.3.43: 


where  N  and  M  are  the  total  degrees  of  freedom  for  displacement  and  pressure 
in  the  entire  structure. 

In  general,  the  integrands  of  (5.3*30)  and  (5-3-38)  are  functions  of 

the  material  coordinates  x3.  It  becomes  very  complicated  when  high  order  of 
shape  function  is  used.  Therefore  the  numerical  integration  technique  must 
be  adopted  in  the  calculation.  However,  when  the  linear  triangular  element 
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is  used,  i^'s,  J,  and  a's  are  constants  within  each  element  at  each 
incremental  intrinsic  time.  This  simplifies  the  calculation  a  great  deal.  On 
the  other  hand  when  a  relatively  fine  mesh  is  used,  the  linear  triangular 
elements  can  provide  the  convergent  results.  The  linear  triangular  elements 
are  still  used  by  many  researchers.  Hence  it  will  be  used  in  this  work.  The 
formulation  of  triangular  element  for  plane  strain  problem  will  be  discussed 
in  detail  in  next  section. 


5.4  Formulation  of  Linear  Triangular  Element  for  Plane  Strain 

The  plane  strain  is  considered  as  a  material  deformation  occurring  in  a 
plane  while  the  deformation  in  direction  normal  to  the  plane  vanl3he3.  Such 
a  deformation  can  be  described  as 


yi  “  >ri(xa*t) 


( i , a  -  1,2)  , 


y,  -  xJ 


(5.4.1) 


The  Cauchy-Green  deformation  tensor  and  its  inverse  and  its  incremental 
form  can  be  calculated  as  follows 

.-1 


[Co8]  " 


and 


CC„83 


C  0 

uv 

0  1 


C  0 

yv 

0  0 


•  CCaB] 


■1 


C  0 

uv 

0  1 


(u, V  -  1,2) 


15.4.2) 


-1 


where  C  ,  C  ,  and  C  are  functions  of  x*  and  x2  only.  The  stresses 
uv  yv  yv 
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related  to  the  x*  direction  from  (5.2.1)  are  the  following 

Ti  j  ■  <Ji  *  “  0  > 

T2 *»  0**-  0  (5.4.3) 

and 

T, o,,-  X  log(J)/J 

where  J  —  1 C  J 1/2  -  |C  i 1/2  remains  the  same  for  both  cases.  Hence  we  need 
1  a8‘  1  uv1 

only  deal  with  the  quantities  related  to  the  xl  and  x*  directions. 

In  the  linear  triangular  element,  the  deformation  can  be  written 
explicitly  as 

-  x1.  V  X1-  »”a*\  (1-1. 2.  a-0, 1,2;  -1,2.  ...6)  , 

y,  -  x*  (u,-0)  (5.4.4) 

where  the  $'s  comprise  a  set  of  constants  related  to  the  nodal  coordinates. 
For  a  typical  triangular  element  shown  in  Fig. 31.  we  denote  the  nodal 
deformation  as 

T 

T 

-  (ultvltu2,va,us,v,)  (5.4.5) 

* 

where  u’s  are  the  displacement  in  xl  direction  and  v's  in  xl  direction, 

T 

respectively.  The  symbol  (•)  denotes  the  transpose  of  a  vector. 
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Then,  we  can  write  the  displacement  in  eq. (5-4.4)  more  explicitly  as 
following 


(5.4.6) 

V  • 


0,-  (x*2)-  x*^))^2^)  (5.4.T) 

G‘“  ^X(3)"  X(2)^2^ 

and  other  a's,  b's  and  c's  can  be  obtained  by  the  subscripts  (1),  (2),  and 
(3)  permuting  in  a  natural  order,  h  is  the  area  of  the  triangle, 


1 

x(1) 

X  * 

x(1) 

A  - 

1 

X(2) 

X  2 

(2) 

(5.4.8) 

1 

x(3) 

X  * 

X(3) 
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a  i 

0 

aa 

0 

a> 

0 

U  -  [1  X1  X1] 

bi 

0 

b* 

0 

b. 

0 

c, 

0 

c2 

0 

Cj 

0 

t 

and 

* 

0 

a, 

0 

a* 

0 

a, 

V  -  [1  X 1  xa] 

0 

b, 

0 

b. 

0 

b. 

0 

Cl 

0 

ca 

0 

c. 

where 

a‘“  (x(2)X 

(3)  " 

X(3)  X 

(2) 

)/(2&) 

The  derivatives  of  displacements  with  respect  to  the  material  system 


Ul,a"  ^iaSn 


(a*1 ,2;  i«1 ,2;  m-1,6) 


(5.4.9) 


where  the  constants  $'s,  in  matrix  form,  can  be  written  as 


and 


r,  1° 


b 

n 

0 


c 

n 

0 


for  m«(2n-l ) 


n  -  1,2  3  (5.4.10) 

for  m  -  2n 


Since  $'s  as  well  as  iji's,  J,  and  a's  are  constants  at  each  intrinsic 
time  step,  the  element  stiffness  (5.3.29)  and  (5.3.30)  become,  in  the  case 
of  linear  triangular  element  as 

Cl ■  a(fj>«  *  4<fjV 

where  the  element  stiffness  matrix  k.  is 

jm 


kjm*  A  UOkk+  25*ak)^k_2^k3),|,il<(Josi)~’,'3i’t,ik(Jok3)] 
and  the  incremental  force  vector 


-  -*ikaQki4 

At  this  moment,  the  degenerate  cases  of  (5.4.10)  and  (5.4.12),  namely 
the  case  under  the  small  deformation,  are  given  bellow.  In  eq.  (5.4.12)  if 
o's  -  0  ,  and  <j>’s  -  $’3,  the  equations  represent  the  small  deformation  of 
plasticity  C55], 
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(5.4.13) 


Dc 


%  “  d(Vex  *  d(fjV  dfj 


(t) 


where 


kj,n* A  ClV*kk*  20(0>*ak)*ik]  • 


“Vp  •  -‘ik^ki1  d* 


(5.4.14) 


Furthermore,  if  f  ■  0,  the  equations  represent  the  small  deformation  of  the 


elasticity  [57], 

CW  <V 


ex 


where  the  element  stiffness  matrix  k.  is 

Jm 

v a  CiOJkk*  20«(ik)*Jik] 

and  (f.)  is  an  external  force  vector. 

J  ex 


(5.4.15) 


(5.4.16) 


5.5  Brief  Discussion  of  the  Calculation  of  Q's 

It  has  been  seen  that  the  effect  of  plasticity  is  represented  by  the 

pseudo-force  df  in  (5.5.12).  In  the  calculation  of  this  force,  dQ  plays  a 
P  QlD 

very  important  role.  Recall  dQ  _  in  (5.2.3) 

ap 


^  fZ  dG(  z-z’  )  dCa6  J_. 

dQa8(z)"  r°  ~ - -  512 


dz  1 


(5.2.3) 


At  each  incremental  step,  dC  .  need3  to  integrate  from  0  to  z,  the  current 

ao 


Intrinsic  time  scale.  To  calculate  dC  .  numerically,  the  equation  for  dC  , 

dp  ci  j 
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at  z.  -  z  ♦  Az  after  knowing  the  value  of  C  at  z.  is  now  derived. 

l  *  i  l  dpi. 

Eq. (5.2.3)  at  z^+1  can  be  written  aa  follows 
,  ,zi»i  aG(z-i')  dCo8  ... 

aWiw«  — - 3?“lz 


Z1  dO(t-z')  fat  dWl-z') 

0  dz  dz'  z ,  dz  dz' 


(5.5.1) 


From  the  mean  value  theorem  and  the  smoothness  of  C  „  ,  the  last  term 

aB 


in  eq. (5.5.1)  can  be  written,  as  a  good  approximation, 

;Zj>A2dG(z-z'  )  dCaS  dCaB  I  ;Zi+AZi  dG(z-z'  )  ^ 

dz  dz'  *  dz  |z-zi  +  1  ^  dz  ^ 


(5.5.2) 


Substitution  of  G(z)  *  T  G  e  arZ  into  (5.5.1)  gives 

r 


a_AZ, 


dC 


aB 


z  +1  Gr(e'arAZi  -  1)}  1-0,1,... 


where  Qa0(O)-O.  (5.5.3) 

Eq. (5.5.3)  tells  that  the  history  dependence  of  the  material  response 
[through  dC  a(z.  ,)]  at  the  intrinsic  time  z.  ,  will  be  determined  by 

Q  a(z  )  plus  the  effect  caused  only  by  the  new  incremental  step  through 
do  i 

dC  - 

-ar  |ztM  ^  tei*i  • 


For  r-1,  i.e.,  G(z)-  G-  e 


-az 


dQ  (2  ) 

aS  i*1 


dC  . 

a  ,  4zi  .  -jfi 

aB  i  dz 


V1 


G,(e'°  Azi  -  1  )  . 


(5.5.4) 
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5.6  The  Iterative  Process  and  Programming  Steps 

After  defining  the  problem,  giving  the  dimension  and  tolerance  of  the 
intrinsic  time,  meshing  the  domain,  and  calculating  the  values  of  <>'3  and 
the  area  of  each  element,  the  incremental  load  step  starts.  For  each 
incremental  step  the  values  of  ' 3  are  updated  and  then  the  stiffness  matrix 

is  calculated.  An  initial  value  Az®  is  assigned  to  the  increment  of 
intrinsic  time  at  every  step.  The  plastic  pseudo-forces  corresponding  to  the 

dz°-are  then  evaluated.  Mow  the  linear  simultaneous  equations  of  incremental 
nodal  displacements  are  set.  After  imposing  the  boundary  conditions,  this 
set  of  equations  is  solved.  Upon  use  of  the  incremental  displacements, 
incremental  strains,  stresses,  and  Cauchy-Green  deformation  tensor  etc.  are 
obtained.  The  incremental  intrinsic  time  is  calculated.  Knowing  the  new 

incremental  intrinsic  time,  the  new  pseudo-force  are  again  obtained.  The 
simultaneous  linear  equations  with  the  3ame  stiffness  matrix  are  solved  and 
new  incremental  nodal  displacements  are  obtained.  The  new  dz  follows  to  the 
new  displacements.  The  iteration  process  is  continued  until  the  difference 
in  dz  between  any  two  consecutive  iterations  is  less  than  some  defined 
tolerance.  Then  the  next  new  incremental  step  is  repeated.  Above  iteration 
procedure  is  described  in  the  flow  chart  as  shown  in  Fig. 32. 

On  the  basis  of  the  formulae  in  section  5.^  and  the  flow  chart  in 
Fig. 32,  a  computer  program  with  Fortran  language  called  FELP  is  developed  to 
analyze  the  plane  strain  of  finite  plastic  deformation  problems.  The 
computer  program  consists  of  five  parts. 


Part  one 


In  this  part  the  dimensioning  of  the  program  variables  is  first  set  up 
and  the  main  control  parameters  (such  as  the  total  number  of  elements, 
nodes,  the  number  of  degree  of  freedom  for  each  node,  etc.)  and  the  material 
constants  (G0,  a.  and  X)  are  input.  The  informations  of  mesh,  namely  the 

coordinates  for  every  node  and  the  nodes  for  each  element  are  then  input. 
The  nodes  of  elements  sure  stored  in  an  array  NQD(N,I)  In  which  N  denotes  the 
number  of  element  and  I  the  number  of  node.  The  array  plays  a  very  important 
role  in  the  connection  of  elements  and  whole  structure.  The  node  and  element 
data  are  defined  by  two  ways  in  the  program,  which  are  automatic  mesh 
forming  in  the  compute  program  and  the  data  inputlng  from  read  statement. 
Former  way  is  efficient  for  relatively  regular  domains  while  the  later  way 
is  for  irregular  ones.  The  necessary  boundary  conditions  are  read  in. 
Finally  all  input  and  calculated  data  are  printed  out  for  checking  and 
recording. 

Part  two 

The  quantities  related  to  the  material  system  are  calculated  in  this 
part,  which  include  the  area  of  the  elements  and  $'s,  i.e.,  a's,  b's,  and 
c,s  defined  in  eq.  (5.^.7).  These  quantities  also  are  primary  variables  to 
the  stiffness  matrix.  To  save  storage  the  stiffness  matrix  is  stored  in  the 
half  bandwidth.  So  the  width  of  the  half  band  is  obtained  in  this  part.  The 
arrays  are  initialized  and  get  ready  for  following  main  calculation. 
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Part  three 

We  start  outer  do-loop  for  the  Incremental  steps  of  load  or 
deformation.  The  values  of  incremental  forces  or  deformations  at  each  step 
are  applied  through  the  boundary  conditions.  For  each  step  the  values  of  f’s 
are  updated.  Then  the  stiffness  matrix  for  each  element  is  obtained  upon  use 
of  eq.  (5.4.1 2)  and  summed  to  the  half  banded  global  stiffness  matrix.  This 
is  accomplished  by  a  subroutine  called  STIFF.  What  follows  next  is  imposing 
the  displacement  boundary  condition  in  the  equation  to  get  the  modified 
global  stiffness  matrix  and  the  modified  force  vector  through  a  subroutine 
BNAR .  The  initial  force  boundary  condition,  i.e.,  node  forces  are  added  to 
force  vector  directly.  It  has  been  seen  from  eq.  (5.4.11)  that  to  complete 
force  vector  needs  to  add  the  pseudo-plastic  force.  We  will  describe  this  in 
inner  do-loop  in  part  5. 

Part  four 

The  inner  do-loop  is  designed  for  the  iteration  process.  At  each 

incremental  step,  an  initial  values  of  intrinsic  time  dz°  is  given.  The 
pseudo-plastic  forces  are  calculated  upon  use  of  the  eq.(5.4.l2)  and  the 
subroutine  PF .  Then  they  are  added  to  the  appropriate  positions  in  total 
force  vector.  Now  the  linear  simultaneous  equations  for  the  nodal 
displacements  are  set.  The  equations  are  solved  by  use  of  the  Gauss- Jordan 
method.  After  obtaining  the  displacement  increments,  new  intrinsic  time  dz 
can  be  obtained  by  calling  a  subroutine  D2.  The  iteration  process  continues 

until  the  value  (dz9-dz)/dz  less  than  the  tolerance. 
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Part  five 


When  inner  do-loop  stops  the  converge  displacement  increments  are 
obtained  at  the  incremental  step.  Then  the  necessary  quantities  such  as 
stains,  stresses,  and  other  related  parameters  through  a  subroutine  called 
RES  can  be  calculated.  If  the  results  for  incremental  values  are  Interested, 
they  are  printed  out.  Otherwise  they  are  added  to  the  total  values  and  the 
next  step  starts.  The  program  consists  of  the  main  program  and  10 
subroutines.  It  will  listed  in  Appendix  A. 

5.7  Numerical  Example 

Previous  analysis  is  now  applied  to  study  a  problem  of  metal  forging  — 
-  upsetting  process  of  a  block  shown  in  Fig. 33.  A  similar  problem  "The 
upsetting  of  a  cylindrical  block  ”  was  taken  up  by  a  Joint  Examination 
Program  of  the  Validity  of  Various  Numerical  Methods  for  the  Analysis  of 
Metal  Forming  Processes  and  discussed  around  table  by  fourteen  groups  in  the 
IUTAM  SYMPOSIUM  TUTZING/ CERMANY  [58].  The  collected  results  showed 
considerable  discrepancies.  However,  what  agreed  in  discussion  was  two 
important  factors  responsible  for  the  discrepancies,  which  are  the 
selections  of  the  elements  and  the  deformation  increments.  What  was 
suggested  was  to  use  finer  element  and  smaller  deformation  Increment  up  to 
0.25S.  To  explore  the  computational  capability  of  the  theory  in  our 
research,  following  data  will  be  taken  in  the  calculation. 

The  dimension  of  original  block  are  set  to  be  2  unit  in  width,  3  in 
height,  and  1  unit  in  third  direction  since  the  plane  stain  is  concerned. 
Sticking,  i.e.,  no  slip  condition,  is  assumed  along  the  tool-work  interface. 
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The  material  constants  G„ ,  X,  and  a  are  set  to  be  1.,  1.,  and  200.  The 

geometrical  dimension  and  the  material  constant  are  taken  in  such  a  way  so 
that  all  quantities  in  the  problem  are  normalized.  G0-X  in  elastic  case 

represents  v-0.25.  Due  to  the  symmetry  of  the  problem,  a  quarter  part  of  the 
block  is  analysed  by  using  the  mesh  division  depicted  in  Fig. 34  in  which  we 
have  143  nodes  and  252  elements.  The  increment  per  each  reduction  (in  hlght) 
step  is  set  to  be  0.001. 

Fig. 35  shows  the  deformation  profiles  for  different  reduction  levels  * 
When  reduction  reaches  0.7,  l.e.,  about  50%,  the  folding  was  observed. 
Fig. 36  shows  the  distorted  grid  vs.  the  original  grid.  The  relatively  rigid 
part  has  been  seen  in  the  up-left  of  the  domain,  i.e.,  the  up-middle  of  the 
block.  Fig. 37  gives  the  bulge  ratio,  w  /w„,  at  the  different  stages  of 

QlcLX 

reduction  in  height.  At  large  reduction,  the  material  gets  softer  in  the 

xl-dlrectlon.  The  variation  of  intrinsic  time  z  respect  to  the  reduction  is 
shown  in  Fig. 38. 

Fig  39  illustrates  the  computed  results  of  upsetting  load  as  a  function 
of  the  reduction  in  flight.  The  variations  of  stresses  with  respect  to  the 
reduction  and  the  intrinsic  time  at  elements  No. 251  and  No. 231  are  shown  in 
Fig. 40-43.  Fig. 40  and  42  give  the  variation  of  stresses  at  element  No. 251  in 

x*-direction  with  respect  to  the  reduction  in  height.  The  results  show  the 
stress  in  this  element  reaches  the  maximum  at  a  critical  value  of  the 
reduction  and  continues  to  decrease  so  that  the  sign  of  the  stress  is 
changed.  However  the  variation  of  stress  at  element  No. 231  which  is  located 
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inner  of  the  domain  does  not  show  the  same  phenomenon.  The  phenomenon 
occurred  at  element  Mo. 251  is  due  to  the  dramatic  distortion  of  the 
structure.  In  Fig.  44-45  we  show  the  distribution  of  stresses  in  and  y2 

direction  along  wi-th  the  width  of  the  upset  block. 

The  convergence  in  this  calculation  in  this  problem  is  excellent.  Under 
the  current  chosen  incremental  reduction,  the  average  number  iteration  for 
each  increment  step  is  two.  This  is  a  beauty  of  the  endochronic  theory  in 
finite  element  method,  because  we  control  dz,  the  Intrinsic  time,  in  the 
iteration  which  plays  the  crucial  rule  in  the  endochronic  theory. 
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Fig.31*  Division  of  mesh 


Fig. 36  Distorted  grid  vs.  U»e 


Fig. 38  Tiie  rariatlon  of  intrinsic  rise 

with  respect  to  reduction  in  tieignt 
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CONCLUSIONS  AND  RECOMMENDATIONS 


Finite  plasticity  has  very  significant  applications  in  engineering 
problems.  With  the  rapid  progress  of  the  modern  technology,  developing  a 
sophisticated  theory  and  applying  it  to  solve  engineering  problem  are  the 
tasks  of  the  science  and  engineer. 

Endochronic  theory  of  plasticity  has  come  to  a  fully  developed  stage 
after  about  15  years'  development  since  1970's.  Its  predictive  power  of  the 
behavior  of  materials,  computational  capability,  application  to  the 
practical  problems  have  been  seen  and  proved  through  analysis  and 
experiments.  In  this  research  we  have  studied  extensively  the  endochronic 
theory  of  plasticity  and  its  practical  applications  under  large 
deformation. Again  we  show  the  contribution  of  endochronic  theory  to  the 
plasticity  regime. 

The  endochronic  theory  of  incompressible  plasticity  was  reviewed  and 
applied  to  the  analysis  of  a  set  of  special  problems  which  have  significance 
in  the  theory  as  well  as  the  practice  of  metal  forming  processes  and 
pressure  vessel  structures.  To  our  knowledge,  this  is  the  first  time  in  the 
field  that  this  set  of  problems  are  solved  in  closed  form  solutions. 

The  theory  was  also  extended  to  model  compressibility  of  plastic 
deformation.  Incompressibility  is  an  idealization  and  simplification  of  the 
reality  and  is  the  key  point  of  the  development  of  the  classical  plasticity. 
We  added  a  term,  representing  the  volume  change,  to  the  free  energy  and 
developed  a  constitutive  equation  of  compressible  plasticity  with  a 
functional  term  which  reflects  the  compressibility  of  materials.  Then  with  a 
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special  form  of  the  function,  we  developed  a  numerical  scheme  with  a  finite 


element  technique  and  a  computer  program  code  for  plane  strain  problems. 
Finally,  an  example  of  metal  forming  of  forging  process  (upsetting  a  block), 
was  analyzed  and  solved  by  the  developed  computer  program.  The  solution 
obtained  by  this  method  is  very  resonable.  The  application  of  endochronic 
theory  to  the  large  plastic  deformation  has  a  great  potential  for 
application  to  more  complicated  engineering  problems. 

The  numerical  algorithm  could  be  extended  to  cover  axial  symmetry  and 
three-dimensional  problems.  Carefully  controlled  experiments  are  needed  for 
further  verification  of  the  validity  of  the  theory. 
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APPENDIX  -  COMPUTER  PROGRAM 


PROGRAM  NAME  FELP 

FUNCTION  ENDOCHRONIC  ANALYSIS  OF  LARGE  PLASTIC 

DEFORMATION  OF  PLANE  STRAIN 

C  DISCRIPTION  OF  IMPORTANT  PARAMETERS  IN  THE  PROGRAM 

C  NEM  -  NUMBER  OF  ELEMENT 

C  NNM  -  NUMBER  OF  NODES 

C  NDF  -  NUMBER  OF  DEGREE  OF  FREEDOM  AT  A  NODE 

C  NPE - NUMBER  OF  NODES  PER  ELEMENT 

C  TH  -  THICKNESS  OF  PLANE 

C  GO , GAMA , CR 1 , AR 1  -  MATERIAL  CONSTANTS  RELATED  TO  KERNEL  FUNCTION 

C  NRMAX.NCMAX  -  MAXIMUM  NUMBER  OF  ROW  AND  COLUM  OF  STIFFNESS  MATRIX 

C  MAXE.MAXD - MAXIMUM  NUMBER  OF  ELEMENT  AND  DISPLACEMENT  B.C. 

C  NHBW - HALF  BAND  WIDTH 

C  ERRO.ERR - TOLERANCE  FOR  INTRINSIC  TIME 

C  NBDF.NBSF  -  NUMBER  OF  NODES  (DISPLACEMENT  AND  FORCE  PRESCRIBED) 

C 

C  DISCRIPTION  OF  IMPORTANT  ARRAIES  IN  THE  PROGRAM 
C  GSTIFT (NRMAX.NCMAX)  BANDED  GLOBEL  STIFFNESS  MATRIX 

C  GF , GFE , GFP(NCMAX )  -  VECTORS  OF  TOTAL,  APPLIED,  AND  PSEUDOFORCE,  RESPECTIVELY 

C  XT, YT(NNM)  —  MATERIAL  COORDINATES  OF  NODES 

C  BI’S,CI’S(NEM)  - <  DIEEERENCE  OF  MATERIAL  COORDINATES  B'S  AND  C'S 

C  PSIT(NEM,2,2,6)  -  VALUES  OF 

C  DXYT.DYXT(NEM)  -  DX/DY  AND  DY/DX 

C  CX,CY,CXY(NEM)  -  CAUCHY -GREEN  DEFORMATION  TENSOR 

C  DCX , DCY , DCXY ( NEM ) - INCREMENTAL  OF  CAUCHY-GREEN  DEFORMATION  TENSOR 

C  TX.TY.TXY(NEM)  -  STRESS 

C  DZO.DZI(NEM)  -  PREVIOUS  AND  PRESENT  INCREMENTAL  INTRINSIC  TIME 

C  BQ’S  Q’S(NEM)  -  QUANTITIES  RELATED  TO  CALCULATION  OF  PSEUDOFORCEO 

C  UT.VT(NNM)  -  DISPLACEMENT  AT  NODES  IN  X  AND  Y  DIRECTION 

C  IBDF ,VBDF(NBDF ) - NODES  AND  VALUE  VECTOR  OF  DISPLACEMENT  B.C. 

C  IBSF.VBSF(NBSF)  -  NODES  AND  VALUE  VECTOR  OF  FORCE  B.C. 

C  NOD(NEM, 3)  -  NODES  OF  ELEMENTS 

C  VOT(NEM)  -  VOLUM  OF  ELEMENTS 

C 

IMPLICIT  REAL*8(A-H ,0-Z) 

DIMENSION  GSTIFT(300,50) ,GF(300) ,GFE(300) ,GFP(300) ,GF1 (300) , 

1  DXY(2 ,2) ,DYX(2 ,2) ,NEDGE( 50) , I  PRINT (30) ,DELT(700) , 

1  BRT(700) ,SEL(700) ,FTU(70Q) ,EDZ( 700) ,  SEL1 (700) , 

1  EDZK260) 

COMMON/UNIT1/B1 1 (260), BI2( 260), BI 3(260), Cl  1 (260), CI2(260), 

1  CI3(260) 

C0MM0N/UNIT2/DCX( 260) ,DCY(260) ,DCXY(260) ,DZ0(260) ,DZ1 (260) , 

1  CX( 260) ,CY( 260) ,CXY( 260) ,TX( 260) ,TXY(260) ,TY( 260) , 

1  DXYT( 260,2,2) ,DYXT(260,2,2) ,QX1 (260) ,QY1 (260) ,  QXY1 (260) , 

1  BQX1 (260) ,BQY1 (260) .BQXY1 (260) ,UT( 150) ,VT(150) , 

1  PSIT(260,2,2,6) 

COMMON/UNIT3/XT(150) ,YT(150),NOD(260,3) ,V0T(260) 
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non 


C0MM0N/UNIT4/IBDF(50) ,VBDF(50) ,IBSF(20) ,VBSF(2Q) 

COMMON/ CONST/NEM,NNM,NDF,NPE 
C 

DATA  TH/1.0D0/ 

DATA  NRMAX.NCMAX.MAXE.MAXD/ 300. 50, 260, 50/ 

DATA  GO , GAMA , CR 1 ,AR1/ 1 .d+0, 1 .0D+0, 1 .0D+0, 200. 0D+0/ 

C 

NDF-2 

NPE-3 

C 

READ(15,*)NNP,ND2 

C  NNP  -  NUMBER  OF  INCREMENTAL  STEPS 

C  NDZ  -  NUMBER  OF  ITERATION  FOR  INTRINSIC  TIME 

C 

READ( 15,*)NPRINT 
DO  9605  I-1.NPRINT 
READ( 1 5 ,* )IPRINT(I ) 

9605  CONTINUE 

C  N PRINT  -  NUMBER  OF  INTERMEDIATE  STEPS  AT  WHICH  RESULTS  PRINTED 

C  IPRINT(NPRINT)  -  STEP  NUMBER  VECTOR 

C 

READ(15,*)IMES 

IMES  -  TYPE  OF  INPUT:  IMES-0  INFORMATION  FOR  NODES  AND  ELEMENTS 

FROM  READ  STATEMENT  IMES-1  FROM  CALCULATION  FOR  REGULAR  DOMAIN 

IF (IMES .EQ . 0)GO  TO  9500 
C 

C  INFORMATION  FOR  RECTAGULAR  BLOCK  IN  FORGE  PROCESS 

C  NX. NX  -  NUMBER  OF  NODES  IN  X  AND  Y  DIRECTION 

C  XL.YL  -  LENGTH  OE  PLATE  IM  X  AND  Y  DIRECTION 

C  I TYPE  -  I TYPE- 1  FIXD  END,  ITYPE-0  SHEAR  FREE  END 

READ( 1 5 ,* )NX,NY .XL, YL .ITYPE 
CALL  NODE(NX,NY,XL,YL,NEQ) 

CALL  DFSF (NX.NY, ITYPE ,NBDF ,NBSF ) 

ND-2*NX+2*NY-4 
DO  260  1-1, NX 
NEDGECl )-I 

NEDGE(NX+NY  +I-2)-NNM-I +1 
260  CONTINUE 

NEDGE(NX*1)-3»NX-1 
NEDGE(2*NX+NY-1 )-(2*NX-1 )*(NY-2)+1 
DO  270  1-1, NY- 3 

NEDG  E  ( NX  >1  + 1 ) -N  EDG  E  ( NX  >1 )  •*■  2  *NX  - 1 
NEDGE(2*NX ♦NY-1 +1 )-NEDGE(2*NX*NY-2+I )-2*NX+1 
270  CONTINUE 
GO  TO  9600 
r 

9500  READ(15,*)NEM,NNM 
DO  9300  I-1.NNM 
READ( 15**)XT(I),YT(I) 

9300  CONTINUE 
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DO  9310  1-1 ,NEM 

READ(15,*)NOD(I,1 ) ,NQD(I,2) ,NOD(I ,3) 

9310  CONTINUE  - 

READ( 1 5, *)NBDF 

DO  9320  1-1  ,NBDF 

READ(15,»)IBDF(I),VBDF(I) 

9320  CONTINUE 

READ(15,»)NBSF 

DO  9400  1-1 ,NBSF 

READ( 1 5 . * )IBSF(I ) , VBSF(I ) 

9400  CONTINUE 
C 

9600  NEQ-NNM*NDF 
WRITE ( 16, 300) 

300  FORMAT (4X, ' ELEMENT’ , 1 0X, ' NODES’ / ) 

WRITE ( 1 6,301 ) (N , (NOD(N ,1 ) ,1-1 ,NPE),N-1 ,NEM) 

301  FORMAT(2X,I5,7X,3I5) 

WRITE( 16,310) 

WRITE ( 1 6 , 305) 

305  FORMAT (4X,’ NODE’ ,15X, ’COORDINATES’/) 

WRITE ( 1 6,306) (N ,XT(N ) , YT(N ) ,N-1 ,NNM) 

306  FORM AT ( 1 X , 1 5 , 1 0X , 2F 1 0 . 4 ) 

WRITE( 16,310) 

310  FORMAT (///) 

WRITE ( 1 6 , 31  5)NBDF 

315  F0RMAT(4X, ’DISPLACEMENT  B.C. ',!$/) 
IF(NBDF.EQ-O)  GO  TO  10 

WRITEC 1 6 , 31 4) 

314  F0RMATC7X,'  POSITION  VALUEV) 

DO  5  I-1.NBDF 
WRITE(16,»)IBDF(I),VBDF(I) 

5  CONTINUE 

10  WRITE( 1 6, 31 0 ) 

WRITE( 1 6, 31 6)NBSF 

316  F0RMAT(4X, 'FORCE  B.C. ’,15/) 

IF(NBSF.EQ.O)  GO  TO  8 
WRITE ( 16,314) 

DO  6  I-1.NBSF 

WRITE ( 1 6 , * )IBSF (I ) , VBSF (I ) 

6  CONTINUE 
C 

8  ERR0-0.01 
WRITE(16,310) 

WRITE ( 16, *) 'GO  -' ,G0 
WRITE( 16,*)' LAMD A  -' ,GAMA 
NHBW-0 

DO  15  N-1.NEM 
DO  15  1-1, NPE 
DO  15  J-1.NPE 

NW-(IABS(NOD(N ,1 )-NOD(N ,  J ) ) *  1 )*NDF 
IF(NHBW.LT.NW)  NHBW-NW 
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15  CONTINUE 
WRITE(16,310) 

WRITE( 1 6,325)  NHBW 
325  FORMAT (MX,' NHBW  -M5) 

C 

DO  20  1-1 ,NEQ 
GFE(I)-0.0 
20  CONTINUE 
C 

DYX(  1 ,1 )-1 . 

DYX(2  ,2)-1 . 

DYX( 1 ,2)-0. 

DYX(2, 1 )-0. 

DXY ( 1 , 1 )  - 1 . 

DXY(2 ,2)-1 . 

DXY (1 , 2)-0. 

DXY(2 , 1 )-0. 

DO  35  1-1 ,NEM 
DZO(I )-0. 

DZ1 (I)-0. 

TX(I )-0. 

TXY(I)-0. 

TY(I)-0. 

BQX1(I)-0. 

BQY1(I)-0. 

BQXY1 (I)-0. 

CX ( I ) - 1 . 

CY  ( I )  —  1 . 

CXY(I)-0. 

DCX(I )-0. 

DCY ( I )-0. 

DCXY(I )-0. 

DO  36  11-1,2 
DO  36  12-1,2 

DXYTd.I1  ,I2)-DXY(I1  ,12) 

DYXTd.I1 ,12)-DYX(I1 ,12) 

36  CONTINUE 
35  CONTINUE 

DO  38  I-1.NNM 
UT(I )-0. 

VT(I)-Q.  ' 

38  CONTINUE 
C 

CALL  BC(TH) 

C 

C**»»»*  STARTING  THE  LOADING  STEPS  ***♦»* 
C 

DO  60  NP-1 ,NNP 
DO  M5  I-1.NEQ 
GF1(I)-0. 

DO  M5  J-1.NHBW 
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GSTIFT(I ,J )-0. 

45  CONTINUE 
IRES-0 
C 

CALL  STIFF ( GST I FT , NRMAX , NCMAX , GO , GAMA ) 

C 

IF  (NBDF.EQ.O)GO  TO  100 
DO  102  1-1 ,NBDF 
IE-IBDF(I) 

VE-VBDF(I) 

CALL  BNDRT (NRMAX, NCMAX, NEQ, NHBW, GSTIFT, GF 1 , IE, VE) 
102  CONTINUE 
100  CONTINUE 
C 

IF(NBSF.EQ.O)GO  TO  105 
DO  25  I-1.NBSF 
II-IBSF (I ) 

GF1 (II)-GFI (II)+VBSF(I ) 

25  CONTINUE 
105  CONTINUE 
C 

C$$$$$  STARTING  THE  INTERATION  FOR  DZ  $$$$$ 

C 

DO  110  NZ-1.NDZ 
DO  115  1-1, NEQ 
GFP(I )-0. 

115  CONTINUE 
C 

CALL  PF(GFP,IBDF , NRMAX, MAXD.CR1 ,AR1 ,NBDF) 

C 

DO  160  N-1 ,NEQ 
GF(N )-GF1 (N)-GFP(N ) 

160  CONTINUE 
C 

CALL  SOL (NRMAX .NCMAX , NEQ , NHBW , GSTIFT , GF , IRES) 
IRES-1 
C 

CALL  DZZ(GF, NRMAX) 

C 

DO  175  N-1 ,NEM 
Z1-DZ1 (N) 

ZO-DZO(N) 

ERR-DABS(Z1-Z0)/Z1 
IF(ERR.GT.ERRO)  GO  TO  1 80 
175  CONTINUE 
ITER-NZ 
C 

GO  TO  200 
180  DO  210  N-1  ,NEM 
DZO(N)-DZHN) 

210  CONTINUE 
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110  CONTINUE 


C 

ENDING  THE  INTERATION  FOR  DZ  $$$$$ 

C 

200  CONTINUE 
C 

CALL  RES ( GO , GF , NRM AX , GAMA ) 

QUANTITIES  RELVENT  TO  BLOCH  FORGING  PROCESS 
DISMAX-XL+UTCNNM) 

DISRT-DISMAX/Xu 
DELT(NP)-VT(1 ) 

BRT(NP)-DISRT 
SEL ( N  P ) -TY ( NEM- 5 ) 

SEL(NP)-TY(NEM-21 ) 

EDZ(NP)-DZ1(NEM-5) 

EDZ1  (NP)-DZI  (NEM-21 ) 

STY-0. 

JJJ-NX-1 
DO  254  I-1.JJJ 
STY-TY(4*I-3)+STY 
254  CONTINUE 

FTU(NP)-STY/JJJ 

PRINTING  INTERM ID I ATE  RESULTS 

DO  2005  I-1.NPRINT 
IF(NP.EQ.IPRINTU))  GO  TO  2006 

2005  CONTINUE 
GO  TO  60 

2006  WRITE(16,310) 

WRITE( 16,340) 

340  F0RMAT(4X, 'ELEMENT  STRAIN’/) 

DO  2000  N-1 ,NEM 
EX-(CX(N)-1 . 0d+0)/2.0d+0 
EY-(CY(N)-1 .0d+0)/2. Od+O 
EXY-CXY(N)/2.0d+0 
WRITE( 1 6,2010)N,EX,EY, EXY 
2000  CONTINUE 
2010  FORMAT ( 2X , 1 5 , 3F  20 . 7 ) 

WRITE(16,310) 

WRITE( 16,600) 

600  F0RMAT(4X, 'ELEMENT  STRESSES  '/) 

DO  220  N-1, NEM 

WRITE ( 1 6,6Q5)N,TX(N),TY(N) ,TXY(N ) 

220  CONTINUE 

605  FORMAT(4X,I5,4F15.7) 

WRITE( 16.J10) 

WRITE( 16,410) 

410  FORMAT (4X, '  NODES  DISPLACEMENTS’/) 

DO  230  N-1 ,NNM 
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WRITE ( 1 6 , 41  5)N  ,UT(N ) ,  VT(N ) 

230  CONTINUE 

415  FORMAT (4X, 1 5,  4X.2F15. 10) 

WRITE( 16,310) 

WRITE ( 1 6 , 420) 

420  FORMAT (4X, *  NODES 
DO  250  I-1.NNM 
NCON-O 
SUMX-O.OD+O 
SUMY-O. OD+O 
SUMXY -O.OD+O 
DO  255  N-1 ,NEM 
DO  255  J-1 ,NPE 
NI--NOD(N ,  J ) 

IF(I.NE.NI)  GO  TO  255 
SUMX  -SUMX  +TX  ( N ) 

SUMY-SUMYTY(N) 

SUMXY -SUMXY  +TXY ( N ) 

NCON-NCON+1 
255  CONTINUE 

SUMX -SUMX/NCON 
SUMY-SUMY/NCON 
SUMXY -SUMXY/NCON 
WRITE ( 1 6 , 605 ) I , SUMX , SUMY , SUMXY 
250  CONTINUE 

WRITE(16,310) 

WRITE ( 16,610) 

610  FORMAT ( 4X, ' EDGE’/ ) 

DO  275  1-1 ,ND 
II-NEDGE(I) 

UD -XT ( I I ) +UT (II) 

VD-YT(II )+VT(II ) 

WRITE(16,*)UD,VD 
275  CONTINUE 
60  CONTINUE 

******  ENDING  LOADING  STEPS  ***** 

PRINTING  QUANTITIES  RELVENT  TO  BLOCH 
WRITE(16,310) 

WRITE(16,610) 

DO  280  1-1 ,ND 
II-NEDGE(I) 

WRITE ( 1 6,*)XT(II),YT(II) 

280  CONTINUE 

WRITE ( 1 6 , 31 0) 

WRITE( 16,620) 

620  F0RMAT(4X, 'RATIO' /) 

DO  290  N-1.NNP 

WRITE ( 1 6 , * )DELT(N ) , BRT(N ) 

290  CONTINUE 
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WRITEC 16,310) 

DO  291  N-1.NNP 

WRITE (1 6,*)DELT(N) ,SEL(N) 

291  CONTINUE 
STOP 
END 
C 

SUBROUTINE  SOL(NRM ,NCM,NEQNS ,NBW , BAND , RHS, IRES) 

C ***SOLVING  A  BAND  SYMMETRIC  SYSTEMS  OF  EQNS 
C ***IN  RESOLVING,  IRES  .GT.  0.  LHS  ELIMINATION  IS  SKIPPED 
IMPLICIT  REAL*3(A-H,0-Z) 

DIMENSION  BAND (NRM ,NCM) ,RHS(NRM) 

MEQNS-NEQNS-1 
IF(IRES.GT.O)  GO  TO  40 
DO  30  NPIV-1.MEQNS 
NPIVOT-NPIV+1 
LSTSUB-NPI V+NBW-1 
IF(LSTSUB-GT.NEQNS)  LSTSUB-NEQNS 
DO  20  NROW-NPI  VOT  .LSTSUB 
C**»INVERT  ROWS  AND  COLUMNS  FOR  ROW  FACTOR 
NCOL-NROW-NPIV+1 

F ACTOR “BAND ( N PI V.NCOL)/ BAND (N PI  V, 1 ) 

DO  10  NCOL-NROW, LSTSUB 

ICOL-NCOL-NROW+1 

JCOL-NCOL-NPIV+1 

10  BAND ( NROW, ICOL) -BAND ( NROW, ICOL) -FACTOR *BAND(N PI V.JCOL) 
20  RHS ( NROW ) “RH S ( NROW )  -F  ACTO  R  *RHS  ( N  PI  V ) 

30  CONTINUE 
GO  TO  70 

40  DO  60  NPIV-1.MEQNS 
NPIVOT-NPIV+1 
LSTSUB-NPI V+N BW-1 
IF(LSTSUB,GT.NEQNS)  LSTSUB-NEQNS 
DO  50  NROW-NPI  VOT,  LSTSUB 
NCOL-NROW^NPI V+1 

FACTOR-BAND (NPI V.NCOL) /BAND (NPIV.1 ) 

50  RHS(NROW)-RHS(NROW)-FACTOR*RHS(NPIV) 

60  CONTINUE 
C***BACK  SUBSTITUTION 
70  DO  90  IJK-2.NEQNS 
NPIV-NEQNS-UJK  +  2 
RHS(NPIV)-RHS(NPIV)/BAND(NPIV, 1 ) 

LSTSUB-NPIV-NBW  +  1 
IF(LSTSUB.LT. 1 )  LSTSUB-1 
NPI  VOT-NPI V-1 
DO  80  JKI -LSTSUB, NPI VOT 
NROW -N  PI  VOT- JKI +LSTSU  B 
NCOL-N PI  V- NROW +  1 
FACTOR-BAND ( NROW, NCOL) 

30  RHS( NROW) -RHS (NROW) -FACTOR *RHS(N PI V) 

90  CONTINUE 
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RHS(1 )-RHS(1 )/BAND( 1 ,1 ) 

RETURN 

END 

C 

SU BROUTINE  BNDR Y ( NRM AX , NCMAX , N EQ , NHBW , S , SL , I E , SVAL ) 
C***THIS  PROGRAM  IMPOSES  THE  PRESCRIBED  B.  C.  ON 
C***THE  SYSTEM  MATRIX (BAND ED  SYMMETRIC  MATRIX) 

C***S  IS  THE  STIFFNESS  MATRIX 
C**»SL  IS  THE  LOAD  VECTOR 

C*»*IE  IS  THE  LABEL  OF  THE  VARIABLE  THAT  IS  PERSCRIBED 
C***SVAL  IS  THE  VALUE  OF  PRESCRIBED  VARIABLE 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  S(NRMAX, NCMAX) ,SL(NRMAX) 

IT-NHBW-1 
I -IE- NHBW 
DO  10  II-1.IT 
1-1*1 

IFQ.LT.1)  GO  TO  10 
J-IE-I  +  1 

SL (I )-SL(I )-S(I ,J )*SVAL 
S(I,J)-0.0 
10  CONTINUE 
S(IE, 1 )  —  1 . 0 
SL(IE)-SVAL 
I-IE 

DO  20  II-2.NHBW 
1-1*1 

IF(I.GT.NEQ)  GO  TO  20 
SL(I )-SL (I )-S(IE,II )*SVAL 
SviE.II )-0.0 
20  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  MV(N.A.MAXE) 

C***THE  INVERSE  OF  MATRIX  A(DIMENSION  IS  N) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(MAXE.MAXE) ,B(200) ,C( 200) 

A( 1 , 1 )-1 .0/A( 1,1) 

NN-N-1 

DO  360  M-1.NN 
K-M+1 

DO  300  1-1 ,M 
C(I)-0.0 
8(1 )-0.0 
DO  300  J-1.M 
B(I)-B(I)*A(I,J)*A(J,K) 

300  CONTINUE 
D-0.0 

DO  310  1-1 ,M 
D-D*A(K,I)*B(I) 
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310  CONTINUE 
D-A(K,K)-D 
A(K,K)-1 .0/D 
DO  320  1-1 ,M 
A(I,K>  —  B(I)/D 
320  CONTINUE 

DO  330  J-1 ,M 
DO  330  1-1 ,M 

330  C(J)-C(J)+A(K,I)«A(I,J) 

DO  340  J-1  ,M 
A(K,  J )— C(  J  )/D 
340  CONTINUE 

DO  350  1-1 ,M 
DO  350  J-1 ,M 

350  A(I,J)-A(I,J)-B(I)»A(K,J) 

360  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  DXDY(N ,DYX,DXY , EUV) 

C*»*TO  CALCULATE  DX/DY  AND  DY/DX 
IMPLICIT  REAL*8(A-H,0-Z) 

C0MM0N/UNIT1 / B1 1 ( 260 ) , BI 2( 260 ) , BI 3 ( 260 ) , Cl  1 ( 260 ) , Cl 2( 260 ) , 
1  CI3(260) 

DIMENSION  BE(3),CE(3),EUV(6),DYX(2,2),DXY(2,2) 

BE(I)-BIKN) 

BE(2)-BI2(N ) 

BE(3)-BI3(N) 

CE(D-CIKN) 

CE(2)-CI2(N) 

CE(3)-CI3(N) 

DO  50  1-1,2 
DYX(I , 1 )-0, 

DYX(I ,2)-0. 

DO  50  J-1 ,3 

IF(I.EQ.I)  K-2»J-1 

IF ( I .EQ.2)  K-2*J 

DYX(I , 1 )-DYX(I , 1 )*BE( J )*EUV(K) 

DYX(I , 2)-DYX(I ,  2)+CE(  J )  *EUV(K) 

50  CONTINUE 

DYX( 1 , 1 )-DYXC 1 , 1  )+1 . 

DYX(2,2)-DYX(2,2)d. 

D-DYX(1 , 1 ) *DYX (2,2) -DYX (1 , 2 ) *D  YX  ( 2 , 1 ) 

DXY(1 ,1 )-DYX(2,2)/D 
DXY(2,2)-DYX(1,1)/D 
DXY(1  ,2)—  DYX(1  ,2)/D 
DXY( 2, 1 )--DYX( 2, 1 )/D 
RETURN 
END 
C 

SUBROUTINE  NODE(NX,NY ,XL, YL,NEQ) 
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C***TO  GENERIZE  TOE  NODES  OF  TOE  ELEMENTS  AND  TOE  COORDINATES  OF  THE  NODES. 
IMPLICIT  R£AL*8(A-H,0-Z) 

COMMON/UNIT1 / B1 1 ( 260 ), BI 2( 260 ),BI 3(260), Cl 1(260) ,CI2(260) , 

1  CI3(260) 

COMMON/UNIT3/XT( 1 50) , YT( 1 50) , NOD ( 260, 3) , V0T(260) 

COMMON/CONST/NEM.NNM.NDF.NPE 

NX1-NX-1 

NY1-NY-1 

NEM-NX1*NY1*4 

NNM-(2*NX-1 )*NY1+NX 

NEQ-NNM*NDF 

NOD( 1 , 1 )-NX+1 

N0D(1 ,2)-2 

N0D(1 ,3)-1 

N0D(2, 1 )-NOD( 1,1) 

N0D(2,2)-N0D(1 ,3) 

NOD(2,3)-2*NX 
N0D(3» 1 )-NOD( 1,1) 

NOD(3,2)-NOD(2,3) 

N0D(3, 3)“2*NX+1 
NOD( 4, 1 )-NOD( 1,1) 

NOD(4,2)-NOD(3,3) 

N0D(4, 3)“N0D( 1,2) 

DO  50  1-1 ,NX-2 
DO  50  J-1 ,4 
DO  50  K-1,3 
II-I*4+J 

NOD(II ,K)-NOD( (1-1 )*4+J,K)*1 
50  CONTINUE 

DO  100  1-1, NY-2 
DO  100  J-1 ,NX1 *4 
DO  100  K-1,3 
II-I*NX1*4-fJ 

NOD(II  ,K)-NOD(  (1-1  )*NX1  *4«-J  ,K)+2*NX-1 
100  CONTINUE 
DX-XL/NX1 
DY-YL/NY1 
DO  150  N-1 ,NY 
DO  150  1-1 ,NX 
NI-(N-1 )*(2*NX-1 )+I 
XT(NI )-(I-1 )*DX 
YT(NI)-YL-(N-1 )»DY 
150  CONTINUE 

DO  160  N-1 ,NY1 
DO  160  1-1 ,NX1 
NI-(N-1  )*(2#NX-1  )+NX+I 
XT(NI)-DX*(2*I-1 )/2 
YT(NI)-YL-DY*(2*N-1 )/2 
160  CONTINUE 
RETURN 
END 
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SUBROUTINE  STIFF (GSTI FT ,NRMAX ,NCMAX , GO .GAMA) 

C#*»TO  FORM  THE  STIFFNESS  MATRIX 
IMPLICIT  REAL«8(A-H,0-Z) 

DIMENSION  GSTIFT (NRMAX , NCMAX ) ,PSI (2  ,2,6) ,BE(3) ,CE(3) . 

1  PSI 1 (2 ,2,6) ,PSI2( 2,2,6) ,ELSTIF(6,6) , DXY( 2,2) ,T(2 ,2) 

1  ,BC(6) 

COMMON/UNIT1/BI 1 (260), BI2( 260), BI 3(260), CI1 (260), CI2( 260), 

1  CI3(260) 

COMMON/ UNI T  2/ DCX ( 260 ) , DCY ( 2  60 ) , DCXY ( 260 ) , DZO ( 260 ) , DZ1 ( 260 ) , 

1  CX(260),CY(260),CXY(260),TX(260),TXY(260),TY(260), 

1  DXYT(260,2,2) ,DYXT(260,2, 2) ,QX1 (260) ,QY1 (260) ,  QXY1 (260) , 

1  BQX1 (260) ,BQY1 (260) .BQXY1 (260) ,UT(150) ,VT(1 50) , 

1  PSI T( 260, 2, 2, 6) 

COMMON/UNIT3/XT(15O),YT(150),N0D(26O,3),VOT(260) 

COMMON/ CONST/NEM , NNM , ND F , N P E 

DO  45  N-1 ,NEM 

VOL-VOT(N) 

BEO)-BII(N) 

BE(2)-BI2(N) 

BE(3)-BI3(N) 

CE(1)-CI1(N) 

CE(2)-CI2(N) 

CE(3)-CI3(N) 

DO  50  1-1 ,2 
DO  50  J-1,2 
DXY(I,J)-DXYT(N,I,J) 

50  CONTINUE 

51- 1. 

52- 0. 

DO  65  IS-1,2 
DO  70  IK-1 ,2 
DO  70  IM-1,3 
M1-2«IM-1 
M2-2*IM 

PSI (IS, IK, Ml )-(DXY( 1 ,IK) *BE(IM)+DXY( 2, IK) *CE(IM) )#S1 
PSI (IS,IK,M2)-(DXY( 1 ,IK) *BE(IM) +DXY( 2, IK) *CE(IM) )*S2 
PSIT(N , IS, IK, Ml )-PSI(IS,IK,M1 ) 

PSIT(N , IS, IK, M2) -PSI (IS, IK, M2) 

70  CONTINUE 

51- 0. 

52- 1 . 

65  CONTINUE 

GJ-DYXT(N ,1,1 )*DYXT(N , 2 ,2)-DYXT(N , 1 , 2 )*DYXT(N ,2,1) 

GG-GAMA*DLOG(GJ)*GO 

T(1 ,1 )-TX(N)*GJ 

T  ( 2 , 2 )  -TY  ( N )  *G  J 

T(1 ,2)-TXY(N)*GJ 

T(2, 1 )-T( 1 ,2) 

DO  80  M-1 ,6 
DO  80  K-1  ,2 
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DO  80  1-1,2 
PSI 1 (K,I ,M)-0. 

PSI2(K,I,M)-0. 

DO  80  IS-1,2 

PSI 1 (K,I ,M)-PSI 1 (K,I ,M)+PSI (K,IS,M) *T(IS,I ) 

1  +PSI (IS,I,M)*T(K,IS) 

PSI 2 (K, I ,M)-PSI2(K,I ,M)+PSI(IS,I ,M)*T(K,IS) 

80  CONTINUE 
DO  85  J-1,6 
DO  85  M-1 ,6 
ELSTIF(J,M)-0. 

DO  90  K-1,2 
DO  90  1-1 ,2 

S-PSKl  ,K,J)*(PSI1  (K,I  ,M)+PSI2(I ,K,M) ) 

SI -PSI (I ,K, J )*(PSI (I ,K,M) +PSI (K,I ,M) ) 
ELSTIF(J,M)-ELSTIF(J,M)-S+GG*S1 
90  CONTINUE 

ELSTIF(J,M)-ELSTIF(J,M)*VOL 
85  CONTINUE 
DO  150  1-1 ,6 
DO  150  J-1 ,6 

ELSTIF(I , J )-ELSTIF(I , J )+(PSI (1,1 ,1 )+PSI (2,2,1))* 

1  (PSI (1,1 , J )+PSI (2 ,2 , J ) )*GAMA*VOL 

150  CONTINUE 
C 

DO  92  1-1 ,NPE 
NR-(NOD(N , I ) — 1 )*NDF 
DO  92  II-1.NDF 
NR  -NR  +  1 

L-(I-1 )*NDF+II 
DO  94  J-1.NPE 
NCL-(N0D(N,J)-1 )*NDF 
DO  95  JJ-1.NDF 
M-( J-1 )*NDF+JJ 
NC-NCL+JJ+1-NR 
IF(NC)  95,95,96 

96  GSTIFT(NR,NC)-GSTIFT(NR ,NC ) +ELSTIF(L ,M) 

95  CONTINUE 
9«  CONTINUE 
92  CONTINUE 
45  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  DZZ(GF.NRMAX) 

C*«*TO  CALCULATE  DZ 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  DYX( 2 ,2) , GF(NRMAX) ,DUV(6) ,PSI (2,2,6), 

1  PSI 1 (2 ,2 ,6 ) , PSI 2(2, 2,6), PSI 3(2, 2, 6), DC (2, 2) 
C0MM0N/UNIT2/DCX( 260 ) , DCY ( 260 ) , DCXY ( 260 ) , DZO( 260 ) , DZ1 ( 260 ) , 
1  CX(260),CY(260),CXY(260),TX(260),TXY(260),TY(260), 


L 
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1  DXYT(260,2,2),DYXT(260,2,2),QX1 (260),QY1 (260),QXY1 (260), 

1  BQX1 (260) ,BQY1 (260) ,BQXY1 (260) ,UT( 150) ,VT( 150) , 

1  PSIT(260, 2,2,6) 

C0MM0N/UNIT3/XK 150) , YT( 150) ,N0D( 260,3) ,V0T(260) 

COMMON /CONST/ N  EM , NNM , ND  F , N  P  E 
DO  150  N-1.NEM 
DO  155  1-1  ,NPE 
NI-NOD(N.I) 

DO  155  J-1  ,NDF 

DUV((I-1 )*2+J)-GF( (NI-1 )*2+J ) 

155  CONTINUE 

DO  160  1-1,2 
DO  160  K-1 ,2 
DYX(I ,K)-DYXT(N ,1 ,K) 

DO  160  M-1 ,6 

PSI (I ,K,M)-PSIT(N ,1 ,K,M) 

160  CONTINUE 
CI-CX(N) 

C2-CY (N ) 

Cl 2-CXY(N ) 

DD-C1 *C2-C1 2*C1 2 
VC1-C2/DD 
VC2-C1/DD 
VC  12— C12/DD 
DO  165  M-1, 6 
DO  166  1-1,2 
DO  166  IA-1,2 
PSI1 (I ,IA,M)-0.0 
DO  166  K-1, 2 

PFT 1 (I ,IA ,M) -PSI 1 (I,IA,M)+PSI(I,K,M) *DYX( K, I A ) 

166  CONTINUE 

DO  167  IA-1,2 
DO  167  IB-1,2 
PSI2(IA,IB,M)«0.0 
DO  167  1-1,2 

PSI2(IA,IB,M)-PSI2(IA,IB,M)+PSI1 (I ,IA ,M) *DYX(I ,IB) 

167  CONTINUE 

DO  168  IA-1,2 
DO  168  IB-1 ,2 

PSI3(IA,IB,M)-PSI2(IA,IB,M)*PSI2(IB,IA,M) 

168  CONTINUE 
165  CONTINUE 

DO  170  IA-1 ,2 
DO  170  IB-1,2 
DC(IA,IB)-0. 

DO  170  M-1 ,6 

DC(IA,IB)-DC(IA,IB)+PSI3(IA,IB,M)»DUV(M) 

170  CONTINUE 

DU  1 -VC  1 *DC (1,1) ♦VC  1 2*DC (1,2) 

DU  2-VC 1 2*DC (1,2) * VC2*DC (2,2) 

DU  1 2-VC 1 *DC (1,2) ♦VC  1 2*DC (2,2) 
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DU  21 —VC  1 2*DC( 1 , 1 )+VC2*DC( 1 ,2) 

DZ-DU1 *DU1 +2. *DU1 2*DU21 +DU2*DU2 
DZ-DSQRT(DZ) 

DZ1(N)-DZ 
OCX ( M ) -DC (1,1) 

DCY (N)-DC(2 ,2) 

DCXY(N)-DC(1 ,2) 

150  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  R ES ( GO , GF , NRM  AX , GAMA ) 

IMPLICIT  REAL*8(A-H ,0-Z) 

C***THE  RESULTS  OF  DISPLACEMENT  AND  STRESSES  FOR  EACH  STEP 

COMMON /UNIT 1/ BI1 ( 260) , BI2<  260 ) , 81 3( 260 ) , Cl 1 (260 ) ,CI2( 260 ) , 

1  CI3(260) 

C0MM0N/UNIT2/DCX( 260) ,  DCY ( 260) ,  DCXY (260) , DZ0(260) ,DZ1 ( 260 ) , 

1  CX(260),CY(260)  ,CXY(260),TX(26Q)  ,TXY(260) ,TY(260) . 

1  DXYT (260,2,2), DYXT ( 260 , 2 , 2 ) , QX 1 (260), QY1 (260),QXY1 (260), 

1  BQX1 (260) ,  BQY1 (260) ,8QXY1 (260) ,UT( 150) ,VT(1 50) , 

1  PSIT(260,2,2,6) 

COMMON/UNIT3/XT( 150) , YT( 150) ,N0D(260, 3) ,V0T(260) 

COMMON/ CONST/NEM , NNM , NDF ,N PE 

DIMENSION  DXY( 2,2) ,DYX(2,2) ,EUV(6) ,GF(NRMAX) ,DUV(6) 

DO  100  N-1.NNM 
UT(N)-UT(N)+GF(2*N-1 ) 

VT(N )-VT(N )+GF(2*N ) 

100  CONTINUE 

DO  110  N-1.NEM 
DO  120  I-1.NPE 
NI-NOD(N.I) 

EUV(2*I-1 )-UT(NI) 

EUV(2*I)-VT(NI) 

DUV( 2*1-1 )-GF(2*NI-1 ) 

DUV(2*I)-GF(NI*2) 

120  CONTINUE 

DO  130  1-1,2 
DO  130  J-1 .2 
DXY(I,J)-DXYT(N,I,J) 

DYX(I,J)-DYXT(N,I,J) 

130  CONTINUE 

GJ-OYX( 1 , 1 )*DYX(2,2)-DYX( 1 ,2)*DYX(2 ,1  ) 

DJJ-O.OD+O 
DO  135  M-1 ,6 

DJJ-DJJ+(PSIT(N , 1 , 1 ,M)*PSIT(N,2,2,M) )*DUV(M) 

135  CONTINUE 

G-GAMA*DLQG(GJ)+GO 

TA-QX1  (N)  ♦G*DCX  ( N )  +GAMA*CX  (N )  *DJJ 

TB-GY1(N)+G*DCY(N)+GAMA*CY(N)»DJJ 

TAB-QXY1 (N)-G*DCXY(N )+GAMA*CXY(N ) *DJJ 

T1J-TX(N)»GJ 
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T2J-TY(N)*GJ 
T12J-TXY(N)*GJ 
D1-DXY(1,1 ) 

D2-DXY( 2,2) 

D12-DXY(1 ,2) 

D21-DXY(2,1  ) 

T1-D 1 *D1 *TA+2.*D1 *D21 *TAB+D21  *D21 *TB 

T2-D 1 2»D 1 2»TA* 2 . *D 1 2*D  2*TAB+D 2*D  2*TB 

T1 2-0 1 2*D 1 *TA+D 1 2*D  21 *TAB+D  2*D 1 *TAB+D 2*D  21 *TB 

F1-0.0D0 

F2-0.0D0 

F12-0.0D0 

F21-0.0D0 

DO  140  M-1,6 

F1-F1+PSIT(N,1 ,1 ,M)*DUV(M) 

F2-F2+PSIT(N,2,2,M)*DUV(M) 

F12-F12+PSIT(N,1 ,2,M)*DUV(M) 

F21-F21+PSIT(N,2,1 ,M)*DUV(M) 

140  CONTINUE 

H1-T1J*F1+T12J*F21 
HI 2-T1J*F12+T1 2J*F2 
H21  -T12J*F1+T2J*F21 
H2«T12J*F12+T2J*F2 
TXJ-T1J+T1-(H1>H1) 

TYJ«T2J+T2-(H2+H2) 

TXYJ-T12J+T12-(H1 2+H21 ) 

CALL  DXDY ( N , D YX , DXY , EU V) 

GJ1-DYXC 1 ,1 )*DYX(2,2)-DYX( 1 ,2)*DYX(2,1 ) 

C 

TX(N)-TXJ/GJ1 
TY(N)-TYJ/GJ1 
TXY(N)-TXYJ/GJ1 
DO  145  1-1,2 
DO  145  J-1,2 
DXYT(N,I,J)-DXY(I,J) 

DYXT(N,I,J)-DYX(I,J) 

145  CONTINUE 

CX(N )-DYX( 1,1) *DYX( 1,1) +DYX( 2,1) *DYX( 2 , 1 ) 

CY(N )»DYX( 1 ,2)*DYX( 1 , 2)+DYX( 2 ,2)#DYX( 2,2) 

CXY(N)-DYX( 1,1) *DYX( 1,2) +DYX (2,1) *DYX( 2,2) 

BQX1 (N )-QX1 (N) 

BQY1 (N)-QY1 (N) 

BQXY1 (N)-QXYI (N) 

110  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  PF(GFP,IBDF,NRMAX,MAXD,CR1 ,AR1 ,NBDF) 

IMPLICIT  REAL»8(A-H,0-Z) 

C***TO  CALCULATE  PLASTIC  FORCES 

COMMON/UNIT2/DCXC260) ,DCY(260),DCXY(260) ,DZ0(260) ,DZ1 (260) , 


1  CX( 260) , CY(260) ,  CXY(260) ,TX(260) ,TXY(260) ,TY(260) , 

1  DXYT(260,2,2),DYXT(260,2,2),QX1(260),QY1(260),QXY1 (260), 

1  BQX1(260),BQY1(260),BQXY1 (260) ,UT( 150) ,VT( 150) , 

1  PSIT(260,2,2,6) 

COMMON/ UNIT 3/XT ( 1 50) , YT( 1 50) ,N0D(260,3) ,V0T(260) 

COMMON/ CONST/ NEM , NNM , NO  F , N  P  E 

DIMENSION  ELQT (2,2) ,ELGFP(6 ) , IBDF(MAXD) .GFP(NRMAX) 

DO  115  N-1.NEM 
VOL-VOT(N) 

ELDZ-DZO(N) 

D1-DXYT(N,1  ,1 ) 

D2-DXYT(N,2,2) 

D12-DXYT(N,1 ,2) 

D21 -DXYT(N ,2,1) 

S-DEXP(-AR1»ELDZ) 

S1-CR1*(S-1  .) 

QX-BQX1 (N)*S+DCX(N)*S1 
QY-BQY1 (N)»S+DCY(N)»S1 
QXY-BQXY1 (N)*S+DCXY(N )*S1 
QX1 (N)-QX 
QY1 (N)-QY 
QXY1 (N)-QXY 

ELQT(1 ,1 )-D1*D1*QX+2.*D21*QXY*DUD21*D21*QY 

ELQT ( 1 , 2 ) -D 1 *D 1 2*QX  +D 1 2«D  21  *QXY  +D 1  *D  2*QXY  *D  2 1  *D  2*Q  Y 

ELQT  ( 2 , 2 )  -D 1 2*D  1 2*QX  *2 .  *D  2*D  1  2*QXY  +D  2»D  2*Q  Y 

ELQT (2,1) -ELQT (1,2) 

DO  120  J-1 ,6 
ELGFP(J)-0. 

DO  120  1-1,2 
DO  120  K-1,2 

ELGFP(J)-ELGFP(J)+PSIT(N,I,K,J)*ELQT(K,I)*VOL 
120  CONTINUE 

DO  125  I-1.NPE 
NI-NOD(N.I) 

DO  130  K-1.NDF 
NII-(NI-1 )*2+K 
DO  135  J-1  ,NBDF 
NJ-IBDF(J) 

IF  (NII.EQ.NJ)  GO  TO  130 
135  CONTINUE 

GFP(NII)-ELGFP((I-1 )*2*K)+GFP(NII ) 

130  CONTINUE 
125  CONTINUE 
115  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  BC(TH) 

C***TO  CALCULATE  BI,CI,AND  AREAS 
IMPLICIT  REAL*8(A-H,0-Z) 

C0MM0N/UNIT1 /B1 1 ( 260), BI2( 260), BI 3(260), Cl  1(260), CI2(  260), 
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1  Cl 3(260) 

C0MM0N/UNIT3/XT(150),rr(150),NOD(260,3),VOT(260) 
C0MM0N/C0NST/NEM , NNM , NDF ,N PE 
DIMENSION  BE(3).CE(3).X(3).Y(3) 

DO  30  N-1 ,NEM 
DO  32  I-1.NPE 
NI-NOD(N.I) 

Xd)-XT(NI) 

Y(I)-YT(NI) 

32  CONTINUE 

AflEA-X(1)*(Y(2)-Y(3))+X(2)»(Y(3)-Y(1))+X(3)*(Y(1)-Y(2)) 
DO  31  I-1.NPE 
J-I  +  1 

IF(J.GT.NPE)J-J-NPE 
K-J  ♦1 

IF(K.GT.NPE)K-K-NPE 
BE(I )-(Y( J)-Y(K) )/AREA 
CE(I)-(X(K)-X(J))/AREA 
31  CONTINUE 
BIKN)-BE(I) 

BI2(N)-BE(2) 

BI3(N)-BE(3) 

CIKN)-CE(I) 

CI2(N)-CE(2) 

CI3(N)-CE(3) 

V0T(N)-AREA/2.*TH 
30  CONTINUE 
RETURN 
END 

SUBROUTINE  DFSF (NX ,NY , ITYPE.NBDF ,NBSF ) 

IMPLICIT  REAL»8(A-H,0-Z) 

COMMON/UNIT4/IBDF(30),VBDF(30),IBSF(2O),VBSF(2O) 

IF(ITYPE.EQ. 0)GO  TO  3000 

NBDF-2»NX+NY-1*NX 

NBSF-0 

DO  3010  1-1, NX 
IBDF(2«I-1 )-I*2-1 
VBDF(2*I-1 )-0.D-00 
IBDF(2*I )-I *2 
VBDF(2»I)  —  1  .0D-03 
3010  CONTINUE 

DO  3020  1-1, NY-1 
IBDF(2*NX+I )-(2*NX-1 )*I *2+1 
VBDF(2*NX+I)-0.D+00 
3020  CONTINUE 

DO  3030  1-1, NX 

IBDF(2*NX+NY-1 >1 )-(2*NX-1 )*(NY-1  )*2+2*I 
VBDF(2*NX-*NY-1  *1  )-0.D*00 
3030  CONTINUE 
GO  TO  4000 
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3000  NBDF-NX+-NY 
NBSF-NX 
DO  4010  1-1 ,NX 
IBSF(I)-2*I 
VBSF(I)-1.0D-03 
4010  CONTINUE 

VBSF(1)-5.0D-04 
VBSF(NX)-5.0D-04 
DO  4020  1-1 ,NY 
I3DF(I)-2*NX»(H1  )  +  1 
VBDFCI )-0.D+00 
4020  CONTINUE 

DO  4030  1-1 ,NX 
I3DF(NY+I)-(NX*(NY-1 )+I)*2 
VBDF(NY+I 5-0.D+00 
4030  CONTINUE 
4000  CONTINUE 
RETURN 
END 


